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ABSTRACT 


An  indirect  adaptive  control  algorithm  has  been  studied 
to  control  physical  systems  with  parameter  uncertainties. 
Although  the  particular  algorithm  investigated  is  applicable 
to  a  wide  class  of  discrete  time  linear  systems,  parameter 
convergence  and  therefore  global  stability  of  the  whole  system 
is  guaranteed  only  if  the  external  excitation  contains  a 
sufficiently  large  number  of  frequency  components. 

This  research  study  has  also  investigated  the  possibility 
of  stopping  the  identification  procedure  when  the  parameter 
error  becomes  sufficiently  small,  so  the  controller  auto- 
matically turns  itself  from  adaptive  to  time  invariant,  while 
still  guaranteeing  global  stability  of  the  closed-loop  system. 
In  this  way  the  adaptive  controller  might  be  activated  only 
when  its  gains  do  not  provide  satisfactory  performances. 

Computer  programs  of  the  indirect  adaptive  control  have 
been  written  for  simulation  purposes  using  both  recursive 
least  squares  and  projection  algorithms. 
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I.   INTRODUCTION 

One  of  the  most  challenging,  interesting  and  active  fields 
of  Automatic  Control  is  Adaptive  Control.   To  implement  high- 
performance  control  systems  when  the  plant  dynamics  are  poorly 
known  or  when  large  and  unpredictable  variations  occur,  the 
control  engineers  prefer  to  use  an  important  class  of  control 
systems  called  Adaptive  Control  systems.   Adaptive  Control 
comes  from  a  desire  and  need  for  improving  performance  of 
complex  engineering  systems  with  large  uncertainties.   It  is 
especially  important  in  systems  with  many  unknown  parameters 
that  are  changing  with  time.   Also,  it  can  be  defined  as  a 
special  type  of  nonlinear  feedback  control,  as  a  nonlinear, 
nonautonomous  dynamic  system. 

An  adaptive  controller  can  change  its  behavior  in  response 
to  changes  in  the  dynamics  of  the  plant  and  the  disturbances. 

The  term  adaptive  control  has  been  used  at  least  from  the 
beginning  of  the  1950 's.   With  recent  advances  in  microproc- 
essor technology,  it  has  become  feasible  to  implement  adaptive 
algorithms  efficiently  in  real  time  at  a  reasonable  cost. 

There  are  three  schemes  for  parameter  adaptive  control : 
gain  scheduling,  model  reference  control  and  self-tuning 
regulators.   The  starting  point  is  an  ordinary  feedback  con- 
trol loop  with  a  process  and  regulator  with  adjustable 
parameters.   But  the  main  problem  is  to  find  a  convenient  way 


of  changing  the  regulator  parameters  in  response  to  changes 
in.  process  and  disturbance  dynamics.   The  following  schemes 
differ  only  in  the  way  the  parameters  of  the  regulator  are 
adjusted. 

A.   GAIN  SCHEDULING 

Gain  scheduling  depends  on  finding  auxiliary  process 
variables  correlated  with  the  changes  in  plant  dynamics.   In 
this  way  it  is  possible  to  reduce  the  effects  of  parameter 
variations  by  changing  the  parameters  of  the  regulator  as 
functions  of  the  auxiliary  variables. 

The  block  diagram  of  this  scheme  is  given  in  Figure  1.1. 
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Figure  1.1.   Block  Diagram  of  Gain-Scheduling  System 


B.   MODEL- REFERENCE  ADAPTIVE  SYSTEMS 

In  this  type  of  adaptive  system,  a  reference  model  speci- 
fies the  desired  performance,  and  tells  how  the  process  output 
should  respond  to  the  command  signal.   A  block  diagram  of 
model  reference  system  is  given  in  Figure  1.2.   As  seen,  from 
this  figure,  the  reference  model  is  part  of  the  control  system, 
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Figure  1.2. 


Block  Diagram  of  the  Model  Reference 
Adaptive  System 


This  adaptive  system  has  two  loops :   the  lower  loop  con- 
tains Regulator  and  Process.   The  upper  loop  adjusts  the 
parameters  of  the  regulator,  in  such  a  way  that  the  error 
between  the  model  output  and  the  process  output  tends  to  zero. 
In  this  type  of  control  system,  the  main  problem  is  to  deter- 
mine the  adjustment  mechanism  so  that  a  stable  system  results, 
which  brings  the  tracking  error  to  zero. 

Model  reference  adaptive  systems  can  be  further  subdivided 
into  two  categories:   direct  and  indirect.   In  indirect  con- 
trol the  plant  parameters  are  estimated  and  the  control  param- 
eters are  adjusted  based  on  these  estimates  so  that  the  overall 
plant  transfer  function  matches  that  of  the  reference  model. 
In  direct  control  no  effort  is  made  to  identify  the  plant 
parameters  but  the  control  parameters  are  directly  adjusted  to 
minimize  the  error  between  plant  and  model  outputs. 

It  turns  out  that  the  model  reference  approach  is  applica- 
ble only  to  plants  with  stable  zeroes.   In  fact  the  only  way 
the  closed  loop  transfer  function  (Plant  &  Regulator)  can 
match  the  one  of  the  model  in  its  poles  and  zeroes,  is  by  re- 
moving the  plant  zeroes  by  cancellation  with  closed  loop  poles. 
This  operation  leads  to  the  presence  of  uncontrollable  or 
unobservable  modes  in  the  closed  loop  systems,  which  can  be 
accepted  only  if  they  are  stable  (i.e.,  their  effect  decays 
to  zero  with  time) ,  and  this  constitutes  the  major  limitation 
for  model  reference  adaptive  systems. 


C.   SELF-TUNING  REGULATORS  AND  POLE  PLACEMENT 

A  third  approach  to  the  adaptive  problem  is  the  self- 
tuning  regulator.   A  block  diagram  of  this  type  control  system 
is  given  in  Figure  1.3. 
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Figure  1.3.   Block  Diagram  of  a  Self-Tuning  Regulator 

The  reference  model  of  the  previous  approach  is  replaced 
by  some  more  general  desired  performances;  such  as  error 
minimization  or  desired  closed  loop  poles. 

There  are  two  loops  in  the  system  configuration.   The 
lower  loop  contains  the  plant  and  linear  feedback  regulator. 
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The  upper  loop  consists  of  a  recursive  parameter  estimator 
and  a  design  calculation  adjusts  the  parameters  of  the 
regulator. 

Also  it  is  possible  to  classify  this  adaptive  control 
scheme  as  direct  and  indirect.   This  classification  depends 
on  the  complexity  of  the  design  calculation  block  that  is 
seen  in  Figure  1.4. 


2   System 


Parameter 
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Design 
Calculation 
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Control 
law 


Figure  1.4 


General  Block  Diagram  of  an  Adaptive 
Control  System 
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The  first  method  is  to  parameterize  the  system  directly 
in  terms  of  compensator  parameters. 

This  research  project  will  concentrate  on  the  last  method, 
that  is,  pole  placement  and  indirect  adaptive  control.   The 
evaluation  of  the  control  law  is  indirectly  determined  on 
the  basis  of  parameter  estimates.   This  method  is  also  called 
explicit,  since  the  design  is  based  on  an  explicit  estimation 
of  the  process  model..  The  parameters  of  the  controller  are 
updated  indirectly  via  estimation  of  the  process  parameters. 
Several  algorithms  to  estimate  the  parameters  of  a  linear 
model  are  available  in  the  literature.   In  this  thesis  two 
well-known  estimation  algorithms  will  be  investigated:   Recur- 
sive Least  Squares  and  Projection.   They  represent  a  tradeoff 
between  complexity  of  computation  and  performances,  in  the 
sense  that  best  performances  (with  relatively  high  complexity 
are  obtained  by  using  Recursive  Least  Squares.   Also,  Block- 
processing  is  investigated  to  determine  the  period  of  the 
adaptation  of  the  control  parameters.   Finite  time  persistency 
of  excitation  is  also  studied  and  simulated  for  indirect 
adaptive  control. 

The  main  reason  of  choice  of  indirect  adaptive  control 
for  this  research  comes  from  the  fact  that  it  is  applicable 
to  nonminimum-phase  systems,  and  there  are  no  limitations  on 
zeroes  of  the  plant.   Therefore,  it  is  more  general  than  model 
reference  adaptive  control . 

Because  of  the  effects  of  the  zeroes  of  sampled  data 
systems  on  adaptive  control  theory,  we  start  investigating 
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the  behavior  of  zeroes  of  sampled  data  systems.   This  thesis 
is  organized  as  follows:   In  Chapter  II  zeroes  of  sampled 
data  systems  are  analyzed  by  stating  results  available  in 
the  literature,  and  simulation  studies.   The  indirect  adaptive 
control  problem  is  presented  in  Chapter  III  and  computer 
simulation  studies  are  given  in  Chapter  IV. 

The  computer  programs  are  presented  in  the  Appendices  B, 
C  and  D.   Also  the  description  of  the  programs  is  given  in 
Appendix  A. 
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II.   ZEROES  OF  SAMPLED  DATA  SYSTEMS 

Many  control  strategies  are  based  on  the  assumption  that 
the  zeroes  of  the  plant  are  on  a  stable  region  of  the  com- 
plex plane  so  that  they  can  be  cancelled  by  precompensation 
or  closed  loop  poles.   One  example  is  the  model  reference 
adaptive  control  mentioned  in  the  Introduction.   However  it 
turns  out  that  in  sampled  data  systems  with  Zero  Order  Hold 
(the  most  popular  mean  of  Digital  to  Analog  conversion)  some 
zeroes  of  the  resulting  discrete  time  plant  are  often  in  the 
unstable  region.   In  this  chapter  we  analyze  the  position  of 
the  zeroes  of  sampled  data  systems,  in  relation  with  the  con- 
tinuous time  plant  dynamics  and  sampling  frequency.   The 
structure,  the  number  of  zeroes  outside  the  unit  disc  and  the 
low  frequency  characteristic  of  the  pulse  transfer  function 
are  examined  from  the  transfer  function  for  an  important  set 
of  continuous  time  processes. 

Finally,  it  is  investigated  that  zeroes  of  the  sampled 
data  systems  are  sensitive  to  high  frequency  poles  present 
in  continuous  time  transfer  function. 

A.   TIME  DOMAIN  ANALYSIS  OF  THE  ZEROES  OF  SAMPLED  DATA 
TRANSFER  FUNCTIONS 

Poles  and  zeroes  are  important  parameters  of  linear  time- 
invariant  systems.   The  zeroes  describe  the  way  the  internal 
variables  are  coupled  to  the  inputs  and  the  outputs.   As 
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described  in  [Ref.  6],  the  unstable  zeroes  limit  the  perfor- 
mance of  control  systems,  since  many  design  techniques  are 
based  on  the  cancellation  of  the  process  zeroes.   This  can  be 
done  provided  the  system's  zeroes  are  stable. 

The  transformation  of  poles  in  the  continuous  time  domain 
to  the  corresponding  discrete  time  is  given  as: 


P  T 
P..  ■+   e  1  (2.1) 


where  T  is  the  sampling  period.   This  transformation  maps 
the  left-half  part  of  the  s-plane  onto  the  unit  disc,  so  that 
stability  is  preserved.   But  there  is  no  simple  transformation 
for.  zeroes  from  continuous  to  discrete  time  domain.   The  type 
of  hold  circuit  affects  the  position  of  the  zeroes.   Most 
digital  control  systems  use  a  zero-order  hold,  and  for  this 
type  of  hold  circuit,  the  effects  are  considered. 

In  the  following  discussion,  the  main  results  are  limit 
theorems,  which  give  the  zero  locations  for  small  and  large 
sampling  periods. 

If  the  continuous  time  transfer  function  G(S)  is  rational 
it  follows  from  Equation  (2.2)  that  H(z)  is  also  a  rational 
function 

H(z)   =   (1-z-1)  l( (ePiT)/(z-ePiT) )ReSpiG(s)/s    (2.2) 


where 
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ReS  .   =   Lim   G(s)/s(s-P.) 
s->Pi  1 


and  P.  are  the  poles  of  G(s)/s. 


The  function  H(z)  has  generically  n-1  zeroes.   For  particu- 
lar values  of  the  sampling  period  some  zeroes  may,  however, 
go  to  infinity,  or  they  may  be  cancelled  by  poles,  i.e., 
hidden  modes.   Hidden  modes  should  not  be  considered  as  zeroes 
of  H(z).   In  fact,  there  are  in  general  no  simple  closed  form 
expressions  for  the  zeroes  of  H.   The  limiting  cases  for  small 
or  large  sampling  periods  can  be  characterized.   That  is  ex- 
plained in  the  following  theorem.   The  major  steps  in  the 
proof  are  given  by  [Ref .  9] . 
Theorem  1; 

Let  G(s)  be  a  rational  function 


(s-z, )  (s-z„)  . . .  (s-z  ) 

G(s)   =   K  -. V? —\ 1 \  (2-3) 

(s-p-^  (s-p2)  • • •  (s~Pn) 


and  H(z)  the  corresponding  pulse  transfer  function.   Assume 

that  m  <  n.   Then  as  the  sampling  period  T  ->  0,  m  zeroes  of 

H(z)  go  to  1  as  exp(z^T)  and  the  remaining  n-m-1  zeroes  of 

H(z)  go  to  the  zeroes  of  B    (z)  where  B  (z)  is  the  polynomial 
J  n-m  n 

defined  as 


B  (z)   =   b,n  n  X  +  b,n  n  2  +  . ..  +  b  n        (2.4) 
n  1  z        2  z  n 


and 
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b,n   =        (-l)k'1£n(^J)  ,    k=l,...,n  (2.5) 

K  1=1  K  * 


The  following  results  can  be  observed  from  the  previous 
theorem: 

1.  The  limiting  zeroes  of  a  pulse  transfer  function  depend 
critically  on  the  pole  excess  of  the  corresponding 
continuous  time  system. 

2.  A  continuous  time  system  with  a  pole  excess  larger 
than  two  will  always  give  a  pulse  transfer  function 
with  zeroes  outside  the  unit  disc  provided  that  the 
sampling  period  is  sufficiently  short.   This  may 
happen  for  quite  reasonable  sampling  periods,  and 
sampled  data  systems  with  unstable  inverses  are  thus 
quite  common. 

An  example  is  given  to  illustrate  the  above  discussion. 

Example  2.1; 

Consider  a  system  with  transfer  function 


G(s)   = 


(s+1)3 


The  corresponding  pulse  transfer  function  can  be  computed  as 


H(z)       = 


b,  z      +   b~z    +   b-. 
(z-e      ) 


where 


b±      =      1    -    (1  +T  +T2/2)e   T 


b2      =       (-2  +T  +T2/2)e    T   +     (2  +T  -T2/2)e    2T 
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2     -2T     -3T 
b3   =   (1  -T  +T  /2)e     -  e 


H(z)  has  a  zero  outside  the  unit  disc  if  0  <  T  <  1.8399. 

The  result  of  Theorem  1  can  be  understood  by  studying  the 
behavior  of  the  continuous  time  transfer  function  G(s)  =  s 
and  its  pulse  transfer  function.   The  reason  why  Theorem  1 
holds  is  because  for  a  sampling  period  T  small  every  system 
behaves  like  G(s)  =  1/s    with  the  number  of  poles  n  and 
number  of  zeroes  m. 

The  pulse  transfer  function  corresponding  to  G(s)  =  s 
is  given  by 

Tn  B  (z) 

H(z)   =   - (2.6) 

n!  (z-l)n 


where  B  (z)  is  given  in  Equation  (2.4). 

The  polynomials  B   are  found  for  some  values  of  n  using 


the  program  given  in  Appendix  A. 


B1(z)   =   1 


B2  (z)   =   z  +1 


2 
B3  (z)   =   z   +4z  + 1 


B4  (z)       =       z3  +llz2   +llz  +1 


Bc(z)       =      z4  +26z3  +66z2  +26z  +1 


5 
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The  polynomials  B   have  zeroes  outside  or  on  the  unit 
circle  for  n  >  2.   The  unstable  zeroes  are  given  in  Table  2.1 


TABLE  2.1 


UNSTABLE  ZEROES  OF  B  (z) 

n 


n        Unstable  Zeroes  of  B  (z) 

—        n- — - 


2  -1 

3  -3.732 

4  -1,  -9.899 

5  -2.322,  -23.20 

Therefore,  it  can  be  noticed  that  there  are  continuous 
time  systems  with  stable  zeroes  such  that  the  corresponding 
pulse  transfer  function  has  unstable  zeroes. 

It  is  possible  to  give  a  complete  characterization  of  the 
zeroes  of  the  pulse  transfer  function  for  small  sampling 
periods.   A  similar  result  for  large  sampling  periods  is 
given  by  Theorem  2 . 
Theorem  2 : 

Let  G(s)  be  a  strictly  proper  rational  transfer  function 
with  G(0)  ^  0  and  ReP .  <  0.   Then  all  zeroes  of  the  pulse 
transfer  function  go  to  zero  as  the  sampling  period  T  goes 
to  infinity. 
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A  lower  limit  on  T  for  stable  zeroes  was  obtained  in 
[Ref.  14],  here  stated  only  for  simple  poles. 


T   >  j   £n[2A(n+l) ]  (2.7) 


where: 


6   =   -  max  ReP.   >   0  (2.8) 


and 


A   =   max|Ai/G(0)  |  (2.9) 

i 

If  G(0)  =  0  it  follows  from  Equation  (2.10); 


P.T 

H(z)   =   G(0)z"1  -Ml-z"1)  I      A.   e  p  T  (2'10) 

i=l         i 
z  - 


The  corresponding  pulse  transfer  function  has  one  zero  at 
z  =  1.   The  behavior  of  the  rest  of  the  zeroes  may  be  more 
complex  as  is  shown  by  Example  2.2. 
Example  2.2: 

Let  continuous  time  transfer  function  be 


G(s) 


[  (s  +  1)2  +1]  (s  +  2) 


Then  the  corresponding  pulse  transfer  function  is 
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—  T  —  T  —  T  — T 

,    .  e       (z-1)  {  (e      +  sinT-cosT)z+e       [1-e       (sin  T+  cos  T)  ]  } 

2z(z-e         )  (z      -2ze        cos  T  +e         ) 


The  zeroes  are  z,  =  1  and 


-2T                 -T 
e    (sin  T+  cosT)  -e 
•?        =   : — . : 

2         -T 

e   +  sin  T-  cos  T 


For  control  system  design  in  discrete  time  domain  it  is 
important  to  know  whether  the  zeroes  of  pulse  transfer  function 
are  inside  the  unit  disc  or  not.   When  all  zeroes  are  inside 
the  unit  disc  the  sampled  system  has  a  stable  inverse  and  all 
its  zeroes  can  be  cancelled.   It  is  therefore  of  interest  to 
find  sufficient  conditions  which  guarantee  that  all  zeroes  of 
a  sampled  transfer  function  are  inside  the  unit  disc.   The 
criteria  are  given  in  Theorem  3  by  [Ref .  6] . 
Theorem  3: 

If  G(s)  is  a  strictly  proper,  rational  transfer  function 
with 

(i)     ReP.   <   0 

(ii)    G(s)   =   0 

(iii)   -7T  <  arg  G(iu))  <  0,   for   0  <  w  <  °° 

Then  all  the  zeroes  of  the  corresponding  pulse  transfer 
function  H(z)  are  stable,  so  that  criteria  for  no  unstable 
zeroes  are  given  both  in  terms  of  G(s)  and  in  terms  of  condi- 
tions on  the  Nyquist  curve  G(ioo). 
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B.   FREQUENCY  DOMAIN  ANALYSIS  OF  THE  ZEROES  OF  SAMPLED 
DATA  SYSTEM 

In  this  section,  the  number  of  zeroes  outside  the  unit 
disc  and  the  low  frequency  characteristic  of  the  pulse  trans- 
fer function  are  analyzed  from  the  transfer  function  of  a 
set  of  continuous  time  processes. 

The  discrete  frequency  response  of  a  continuous  process 
G*(co)  can  be  derived  from  Equation  (2.11)  by  substituting 
z  =  exp(jooT)  : 


(z-o-,  )  ...  (z-a   -,  ) 

G(z)   =   K  -t —. -, ^i-  (2.11) 

(z-p1)  . . .  (z~Pn) 


or  it  can  be  expressed  by  the  holding  device  transfer  function 
H(co)  and  process  transfer  function  G(cu).   Then  it  becomes: 


1 
G*(u>)   =  ±       I       H(w+kft)  G(w+kft)  (2.12) 

k=-°° 


where  k  is  an  integer  and  0.   =  2tt/T.   G(lo)  is  the  frequency 
response  of  the  continuous  time  system  G(s). 

For  zero-order  hold,  the  holding  device  frequency 
response  is 


i    -JwT 

H(tu+kQ)  =   .  ,  e^ .  n .  (2.13) 

j  (to  +k^) 


Substituting  Equation  (2.13)  into  Equation  (2.12)  and  taking 
into  account  that  G(-co)  =  G(oj),  where  G  is  the  conjugated 
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value   of   G,    we   obtain 


G*(co)       = 


1  -e 


-jcoT 


jcoT 


G(co)  [1  -f  (co)    +£  (co)  ] 


(2.14) 


where : 


f  (00) 


co       G(ft-co) 


J2  -co      G  (co) 


(2.15) 


£  (co) 


CO 


I  I 


GirQ.  +  co)       G(  (r+1)  fi  -co) 


G(co)       ^-,         rft  +co  [r+I)fl  -co 


}  ~ 0   (2.16) 


and  r  is  an  integer. 

The  number  of  the  zeroes  outside  the  unit  disc  are  given 
by  Theorem  4  [Ref .  10] . 
Theorem  4 : 

If  a  continuous  process  satisfies  Equation  (2.16)  and  the 
condition 


G(fl-co) 
G(co) 


<   1 


(2.17) 


the  number  of  poles  of  G(z)  outside  the  unit  disc  is  zero  and 
the  phase  angle  <p(ft/2)    =   arc  G(Q/2)    is  between  0  and  -tt  then 
its  pulse  transfer  function  possesses  no  zeroes  outside  the 
unit  disc.   If  the  value  of  phase  angle  is  between  -tt  and  -2tt 
one  of  the  zeroes  lies  outside  the  unit  disc,  etc. 
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(z-a  )  ...  (z-a. ) 
G*(z)   =   K*  — 


-T/T,  T/T, 

(Z-e      )  ...  (z-e    b) 


Thus,  inverse  stable  continuous  processes  frequently 
possess  inverse  unstable  pulse  transfer  functions.   For 
instance  if  in  a  transfer  function  satisfying  the  conditions 
given  by  Equations  (2.16)  and  (2.17),  t.  >  T  and  T.  >  T  are 
real  and  positive  for.  every  i  and  the  orders  of  its  denominator 
and  numerator  differ  by  more  than  two,  then  usually  the  pulse 
transfer  function  is  inversely  unstable. 

Also,  from  Theorem  4,  the  phase  of  the  continuous  frequency 
response  at  the  Nyquist  frequency  fi/2  is  directly  related  to 
the  number  of  unstable  discrete  time  zeroes.   High  frequency 
dynamics  (where  high  frequency  is  intended  as  above  the  Nyquist 
frequency)  influence  the  phase  at  the  Nyquist  frequency,  so 
that  it  can  cause  the  number  of  unstable  discrete  time  zeroes 
to  increase.   This  high  frequency  pole  may  come  from  the  struc- 
ture of  the  measuring  devices,  actuators  or  other  hardware  parts 
of  the  physical  realization.   Usually,  we  neglect  these  terms 
in  the  transfer  function.   The  following  example  can  give  an 
idea  about  this  problem. 
Example  2.3: 

Let  the  transfer  function  of  the  plant  be 


H(s)   = 


2   i 
s  - 1 


which  has  a  pulse  transfer  function  given  by 
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„,  .      .54308z  +  .54307 
H(z)   =   — ~ 

z  -3.08616Z  +1 


with  one  zero  at  z  =  -0.99  8. 

By  adding  the  extra  dynamics  D(s)  as 

2 
w 

D(s)   = 


2  2 

s   +2^w  +w 
n    n 


the  zeroes  of  the  pulse  transfer  function  of  the  entire  sys- 
tem for  c  =  0.9  and  some  w   values  are  given  in  Table  2.2 
and  plotted  in  Figure  2.1. 


TABLE  2.2 

ZEROES  LOCATIONS  FOR  DIFFERENT  w   VALUES 

n 

Zeroes 

-0.01349725  -0.1775817 

-0.009074569  -0.1386477 

-0.005823303  -0.1097788 

-0.003552481  -0.08844686 

-0.002030142  -0.07263207 

-0.001109288  -0.06071956 

-0.0004927621  -0.05170227 

-0.0002742233  -0.04454243 

-0.00009020962  -0.03895098 

-0.00006002390  -0.03435726 

-0.00004560289  -0.03068041 

-0.00003143626  -0.02745582 

-0.00005455861  -0.02487628 

-0.0001198984  -0.02266212 
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w 

n 

7 

-3.162004 

8 

-2.866706 

9 

-2.633762 

10 

-2.447441 

11 

-2.296317 

12 

-2.172057 

13 

-2.068539 

14 

-1.981240 

15 

-1.906793 

16 

-1.842641 

17 

-1.786843 

18 

-1.737905 

19 

-1.694647 

20 

-1.656146 
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We  notice  that  as  w   decreases,  the  unstable  zero  gets 
more  and  more  unstable. 
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III.   INDIRECT  ADAPTIVE  CONTROL  FOR  DISCRETE  TIME  SYSTEMS 

As  mentioned  in  the  Introduction,  the  main  advantage  of 
indirect  adaptive  control  systems  comes  from  the  fact  that  it 
is  applicable  to  discrete  time  systems  with  unstable  zeroes. 
It  is  assumed  that  the  system  is  as  given  in  Figure  3.1, 


Figure  3.1.   Single  Input — Single  Output  Plant  System 


The  pulse  transfer  function  is  given  by 


H(z) 


r(z) 
p(z) 


(3.1) 


r(z)  and  p(z)  being  polynomials  given  by 


p(z)   =   zn  +  p^11"1  +  ...  +  pn  (3.2) 


-I  s* 

r(z)   =   r,z     +  r2z  "   +  ...  +  rn      (3.3) 
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The  degree  of  the  denominator  polynomial  p(z)  determines  the 
system  order  n.   Also,  if  we  assume  the  plant  to  be  causal, 
the  polynomial  r(z)  has  degree  at  most  n-1.   The  sequences 
u(t)  and  y(t)  denote  the  system  input  and  output,  respectively 
An  alternative  way  of  describing  the  plant  is  by  its  differ- 
ence equation  written  in  operator  form,  as 

p(D)y(t)   =   r(D)u(t)  (3.4) 

where : 


p(D)   =   1  +  p  D  +  ...  +  pRDn  (3.5) 


r(D)   =   r-^  +  ...  +  rRDn  (3.6) 


and  D  =  q    is  the  backward-shift  operator,  defined  as 
Dy(t)  =  y(t-l) . 

In  the  adaptive  control  problem  the  parameters  in  p(D)  and 
r(D)  are  assumed  to  be  unknown.   Only  the  system  order  n  is 
assumed  to  be  known  to  the  designer.   Hence,  two  problems 
appear  at  this  point: 

1.  Parameter  estimation;  and 

2.  Compensator  structure. 

The  general  block  diagram  of  the  indirect  adaptive  control 
problem  is  given  in  Figure  3.2.   Estimation  of  the  parameters 
of  the  plant  is  mentioned  in  Chapter  III.C.   The  compensator 
structure  is  derived  from  the  pole  placement  problem  for 
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Figure  3.2.   General  Block  Diagram  of  Indirect  Adaptive 
Control 


linear  systems;  and  it  is  computed  on  the  basis  of  the  esti- 
mated plant  parameters . 

In  the  next  section,  we  discuss  the  pole-placement 
problem  which  will  be  the  basis  for  the  determination  of  a 
suitable  controller  structure. 

A.   POLE  PLACEMENT  DESIGN  BASED  ON  INPUT-OUTPUT  MODELS 

The  purpose  of  pole-placement  by  state-feedback  is  to 
determine  a  feedback  controller  so  that  all  poles  of  the 
closed-loop  system  assume  prescribed  values.   This  can  be 
easily  achieved  if  all  state  variables  of  the  plant  are 
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measured.   In  general,  the  states  may  not  be  directly  avail- 
able.  However,  under  certain  observability  conditions  the 
state  can  be  observed  on  the  basis  of  input-output  measure- 
ments.  For  this  reason  observers  are  used  to  estimate  the 
state  of  any  observable  realization.   Along  these  lines,  the 
controller  can  be  considered  as  composed  of  an  observer  and 
a  gain  matrix. 

Consider  a  state-space  representation  of  the  plant  assumed 
to  be  controllable  and  observable: 

x(t+l)   =   $x(t)  +  Tu(t)  (3.7) 

y(t)   =   Cx(t)  (3.8) 

with  <£> ,  r,  C  matrices  of  dimensions  n  xn,  n  xl  and  1  xn, 
respectively. 

The  block  diagram  of  the  controlled  system  combined  with 
the  observer-controller  is  given  in  Figure  3.3.   V(t)  is  an 
external  input  which  will  be  discussed  in  Chapter  IV.   If 
we  restrict  ourselves  to  the  class  of  linear  control  inputs 
we  can  write: 

u(t)   =   -  Lx(t)  +  v(t)  (3.9) 

as  discussed  in  [Ref .  3] ,  where  L  is  a  constant  matrix  as 


L   =   [L  1.  ...  1  ] 
12       n 
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v(t) 


+ 


u(t) 


6 


PLANT 


STATE 
ESTIMATOR 


y(t) 


Figure  3.3.   Block  Diagram  of  the  Controlled  System 

and  x(t)  indicates  the  estimated  values  of  the  actual  states 
of  the  plant.   It  is  a  well-known  result  in  systems  theory 
[Ref.  3]  that  we  may  define  the  observer  dynamics  as 


x(t+l)   =   <J>x(t)  +  ru(t)  +  K[y(t)  -Cx(t)] 


(3.10) 


where  K  is  a  matrix  such  that 


K1   =   [kx   k2 


k  ] 

n 


and  $-KC  is  a  matrix  with  eigenvalues  inside  the  unit  circle 

Rearranging  Equation  (3.10)  we  obtain  the  state-space 
dynamical  equations  of  the  controller  as  a  two  input — one 
output  linear  system. 


x(t+i: 


C$-KC]x(t)  +  ru(t)  +  Ky(t) 


(3.11) 
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u(t)   =   -  Lx(t)  +  v(t)  (3.12) 

To  convert  the  state-space  representation  of  the  pole- 
placement  problem  to  a  transfer  function  (using  polynomials) 
form,  taking  the  z  transform  of  Equations  (3.11)  and  (3.12) 
we  obtain 


X(z)   =   (zI-Q)  1TU(z)  +  (zI-Q)  1KY(z)  (3.13) 


U(z)   =   -L(zI-Q)-1rU(z)  -L(zI-Q)~1KY(z)  +V(z)      (3.14) 


where  we  define  the  matrix  Q  =  $-KC.   Also  we  can  define  the 
rational  functions 


-Lizl-O)'1?      -      L   adj  (2I~Q)T 
L(zl  Q)   r   -     det(zI-Q) 


k(z) 
q(Z) 


(3.15) 


and 


L(ZI-Q)_1K   =   "La/!iZ''°!K 

det(zl-Q) 


h(z) 


q(z) 


(3.16) 


where  q(z),  the  observer  polynomial,  is  an  arbitrary  stable 
polynomial.   Arbitrariness  of  q(z)  is  guaranteed  by  the 
assumption  of  the  plant  being  observable,  while  h(z)  and 
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k(z),  the  controller  polynomials,  have  to  be  computed,  on  the 
basis  of  the  plant  dynamics  and  q(z). 

Thus,  the  structure  of  the  control  input  u  can  be  described 
as  : 


u(2)    =   iHff  u(z)  +Htf  Y(z)  +  v(z)  (3-17) 


or  it  may  be  expressed  in  a  difference  operator  form  as 

q(D)u(t)   =   K(D)u(t)  +  h(D)y(t)  +  q(D)v(t)         (3.18) 

To  make  the  notation  more  attractive,  let  us  write 
Equation  (3.18)  as 


U(t)   =   qW  U(t)  +  qW  Y(t)  +  V(t)  (3'19) 


The  block  diagram  of  the  above  controlled  system  is  in  Figure 
3.4.   Hence,  the  closed-loop  system  defined  from  external 
input  v(t)  to  the  output  y(t)  has  transfer  function 


y(t)   "   p(D)[q(S)-k(D)]  -r(D)h(D)  V(t)  (3-20) 


In  the  pole-placement  problem,  the  goal  is  to  determine 
the  compensator  parameters  (k,h,q)  so  that  the  closed-loop 
poles  are  as  we  desire, 


y{t)    =    ^ihyv{t)  (3-21) 
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v(t) 


Figure  3.4.   Transfer  Function  Form  of  the  Closed- 
Loop  System 


p*(D)  being  an  arbitrary  stable  polynomial  in  the  difference 
operator  D.   Roots  of  p*(D)  are  the  new  assigned  poles  of  the 
closed- loop  system. 

By  equating  Equations  (3.20)  and  (3.21),  we  obtain  that 
k,  h,  q  have  to  satisfy  the  polynomial  equation 


p(D)[q(D)  -k(D)]  -  r(D)h(D)   =   q(D)p*(D) 


(3.22) 


which  can  be  written  as 


k(D)pCD)  +  h(D)r  (D)   =   F(D) 


(3.23) 


where : 


F(D)   =   q(D) [p(D)  -  p*(D)] 
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Therefore,  the  problem  of  determining  a  suitable  compen- 
sator for  pole-placement  is  equivalent  to  solving  the  poly- 
nomial equation  (3.23)  which  is  known  by  the  name  of 
Diophantine  equation.   The  conditions  by  which  this  equation 
can  be  solved,  together  with  the  method  of  solution  itself, 
are  given  in  the  following  section. 

B.   DIOPHANTINE  EQUATION 

The  general  form  of  the  Diophantine  equation  in  the 
unknown  polynomials  K(z)  and  H(z)  is  given  by: 

F(z)   =   k(z)P(z)  +  h(z)R(z)  (3.24) 


with 


k(z)   =   k   +  k,z  +  ...  +  k  zm,   k  ±    0       (3.25) 
o     1  mm 


h(z)      h   +  h, z  +  ...  +  h  z  (3.26) 

o     1  m 


P(z)   =   P   +  p, z  +  ...  +  p  zn  (3.27) 

o    *±  ^n 


R(z)   =   r   +rnz+...+rzn  (3.28) 

o     1  n 


and 


where : 


F(»   =   fQ  +  f,z  +  f2z2  +  ...  +  fn+mzn+m     (3.29) 
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k.,  h.,  p.,  r.  and  f.  are  constant,  not  necessarily 
11*11       1  -1 

all  nonzero. 

In  the  general  pole  placement  setting,  the  polynomials 
p(z),  r(z),  F(z)  are  given,  while  h(z),  k(z)  are  unknown. 

In  the  rest  of  this  section,  we  determine  the  conditions 
by  which  the  Diophantine  equation  can  be  solved  and  the  method 
of  solution. 

By  substituting  Equations  (3 . 25) - ( 3 . 29 )  into  (3.24),  we 
obtain : 


f+f ,  z  + . .  .   +f    ,    zn+m      =       (p     +p,  z+...+pzn)(k     +k,z   +    ... 
o        1  n+m  ^o      ^1  *n  .  o        1 

+   k    zm)    +    (r     +  r,z  +  .  .  .   +r    zn)  (h      + 
m  o         1  no 

h,  z    +    ...    h    zm)  (3.30) 

1  m 


and  equating  the  coefficients  of  the  same  power  of  z  yields 
the  linear  relation: 


CTS    =   [f    f ,   f 0  . . .   f  ^  ]  (3.31) 

m       o    1    2       n+m 


where  the  vector  C  is  defined  as: 


C   =   [k    h    k,   h,   .  .  .   k    h  ]T  (3.32) 

o    o    1    1        mm 


and  the  matrix  S   as: 

m 
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m 


Po   ?! 


ro   rl- 


0    P. 


o 


Pn-1   pn 


0     0.  .  .0 
0     0...0 


0    0 


r   ,   r      0    0 
n-1    n 


.p   0   p   ,   p    0 
^n-2   ^n-1   ^n 


r   0   r   -.   r    0 
n-2    n-1    n 


n 


n 


(3.33) 


Equation  (3.31)  is  a  linear  algebraic  equation  in  the  unknown 

vector  C.   The  matrix  S   consists  of  m+1  block  rows;  each 

m 

block  row  has  two  rows  and  can  be  obtained  by  shifting  its 
previous  block  row  to  the  right  by  one  column.   It  is  a 
2 (m+1)  x  (n+m+1)  matrix. 

In  order  for  a  solution  to  exist,  the  matrix  S   has  to 
be  full  column  rank  [Ref .  5] .   This  can  be  satisfied  if 
2  (m+1)  21  n+m+1,  or  m  >_  n-1. 

When  S   is  a  square  matrix  (m  =  n-1) ,  it  is  called  the 
Sylvester  matrix  of  P  and  R,  which  has  nonzero  determinant 
provided  polynomials  P  and  R  are  mutually  coprime,  i.e.,  they 
do  not  have  common  factors  [Ref.  6] . 

We  can  summarize  the  steps  to  find  h  and  k,  mentioned 
above,  as: 

1.  Let  n  be  the  order  of  the  system. 

2.  Choose  q(z)  to  be  an  arbitrary  polynomial  of  degree  n. 
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3.  Form  the  matrix  S   using  P(z)  and  R(z) . 

4.  Set  up  the  polynomial  F(z)  as: 

F(z)   =   q(z)  [p(z)  -  p*  (z) ]  . 

5.  Solve  for  parameters  of  h(z)  and  k(z). 

6.  Set  up  the  polynomials  h(z)  and  k(z),  or  h(D)  and 
k(D)  using  the  relation  D  =  z~l. 

C.   PARAMETER  ESTIMATION 

As  mentioned  above,  in  indirect  adaptive  control,  we 
have  to  estimate  the  parameters  of  the  plant,  where  the 
parameters  are  the  coefficients  of  the  transfer  function 


H(z)   =   £-^v  (3.34) 

p(z) 


We  can  compute  the  controller  parameters  h(D)  and  k(D) 
from  the- Diophantine  equation  on  the  basis  of  the  estimated 
plant  (say  r ,p) . 

A  large  number  of  different  identification  methods  are 
available.   In  the  literature,  one  broad  distinction  is 
between  on-line  methods  and  off-line  methods.   In  the  off- 
line case,  it  is  presumed  that  all  data  are  available  prior 
to  analysis.   Consequently,  the  data  may  be  treated  as  a 
complete  block  of  information,  with  no  strict  time  limit  on 
the  process  of  analysis.   In  contrast  to  the  off-line  case, 
the  on-line  case  deals  with  sequential  data,  which  requires 
the  paraemter  estimates  to  be  recursively  updated  within  the 
time  limit  imposed  by  the  sampling  period. 
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As  mentioned,  on-line  estimation  schemes  produce  an 
updated  parameter  estimate  within  the  time  span  between 
successive  samples. 

Also,  the  on-line  methods  are  the  only  alternative  if 
the  estimation  is  going  to  be  used  in  an  adaptive  controller 
or  if  the  process  is  time-varying.   In  this  thesis,  we  con- 
sider two  classes  of  estimation  algorithms: 

1.  projection; 

2.  recursive  least-squares. 

Before  proceeding,  it  can  be  said  that  input-output 
characteristics  of  a  wide-class  of  linear  and  nonlinear 
deterministic  dynamical  systems  can  be  described  by  a  model 
expressed  in  the  following  form: 


y(t)   =   <Mt-l)Teo  (3.35) 


where  y(t)  denotes  the  system  output  at  time  t,  $(t-l) 
denotes  a  vector  given  as : 

<Mt-l)T   =   [y(t-l)  y(t-2)  ...  y(t-n)  u(t-l)  ...  u(t-n)] 

(3.36) 

and  G  denotes  the  parameters  of  the  plant,  described  as: 


T 

)Q   =   {-Pl,-p2, . . . ,pn,r1,r2, . . . ,rm)  (3.37) 


The  following  example  illustrates  the  representation  of 
the  plant  as  given  in  Equation  (3.35)  . 
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Example  3.1: 

Suppose  the  pulse  transfer  function  is 


H(z)   = 


z  +1 


z   +2z  +3 


It  is  also  expressed  in  the  following  form; 


HCz"1)    '=      — Z        +7 


-1  -2 

1  +2z    x  +3z 


H(q    1) 


-1  _,_    -2 

__2 ±2. 

-1  -2 

1   +2q    X   +3q 


Hence,  the  difference  equation  becomes 


y(t)   =   -2y(t-l)  -  3y(t-2)  +  u(t-l)  +  u(t-2) 


which  can  be  expressed  as: 


y(t)   =   [y(t-l)  y(t-2)  u(t-l)  u(t-2)] 


-2 

-3 

1 

1 


The  rest  of  this  section  illustrates  the  estimation 
methods.   The  projection  algorithm  and  the  recursive  least-squares 
algorithm  are  analyzed  sequentially. 
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( I )   Projection  Algorithm 

By  the  Projection  Algorithm,  the  sequence  of  estimates 
(t)  is  recursively  computed  as: 


9(t)   =   0(t-l)  +  a°(tT1) [y(t)-$(t-l)T6(t-l)]    (3.38 

c  +<D(t-l)  <S>(t-l) 


with  arbitrary  initial  estimate  9(0)  and  c>0;0<a<2. 
The  detailed  derivation  steps  are  given  by  [Ref .  2] . 

This  algorithm  is  also  known  as  the  normalized  least-mean- 
squares  (NLMS)  algorithm.   The  algorithm  results  from  the 
following  optimization  problem:   Given  6(t-l)  and  y(t0, 
determine  9(t)  so  that 

2 

J   =  j\   |9(t)  -9(t-l) | |  (3.39) 

is  minimized  subject  to 


y(t)   =   $(t-l)T9(t)  (3.40) 


The  main  properties  of  the  projection  algorithm  are  the 

following : 

1.   Estimation  of  9  at  time  t,  9 (t)  is  always  closer  to 
the  actual  value  of  9  than  the  preceding  estimated 
value  9 (t-1)  . 


(t)  -e0||  <  ||9(t-i)  -e0||     I | e  c  o )  -e0||;   (3.4i 


t  >  i 
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2.   When  the  time  goes  to  infinity,  the  error  between  the 
actual  output  of  the  system  and  its  predicted  value 
will  converge  to  zero. 


lim  Sit) =   0  (3>42) 

t->~   [C  +$  (t-1)  ^  (t-1)  ]   7 


where 


e(t)   =   y(t)  -  $(t-l)T9(t-l) 


3.   lim   |6(t)  -9(t-k)|    =   0   for  any  finite  k       (3.43) 

t->oo 


The  adaptation  rate  converges  to  zero  as  t  -*■  °°. 

(II)   Recursive  Least  Squares  Algorithm 

According  to  Gauss  the  principle  on  which  the  least  square 
estimates  are  based  is  that  the  unknown  parameters  of  the 
process  should  be  chosen  in  such  a  way  that  the  sum  of  the 
squares  of  the  differences  between  the  actually  observed  and 
computed  values  multiplied  by  numbers  that  measure  the  degree 
of  precision  is  a  minimum. 

The  recursive  least  squares  algorithm  is  given  by  the 
following  equation  in  [Ref.  2]. 


(t)   =   6(t-l)  +  P(t-2H(t-l)[y(t)-*(t-l)T9(t-l)]        (3>44 

1  +$  (t-1)  P(t-2)<Mt-l) 


for  t  >  1  and 
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P(t-l)   =   p(t-2)  P(t-2)Mt-l)Mt-l)TP(t-2)  (3>45) 

1  +$  (t-1)  P(t-2)$(t-l) 


with  0(0)  given  and  P(-l)  is  any  positive  definite  matrix 

P  . 
o 

The  algorithm  results  from  the  minimization  of  the 
following  quadratic  cost  function: 


JN(6)   =  \      I    {y(t)-0(t-l)T9}2  +  j(e-9(0)T  Po"1(9-6(0)) 

(3.46) 
We  can  observe  that  the  cost  function  represents  the  sum  of 
squares  of  the  output  prediction  error  e(t). 


e(t)   =   y(t)  -  $(t-l)T6  (3.47) 


The  second  term  of  the  right-hand  side  of  the  cost  function 

accounts  for  the  initial  parameter  estimates  weighted  by  the 

matrix  P  .   In  this  way.  we  can  consider  P   (the  initial  con- 
o  J  o 

dition  of  P(t)  in  Equation  (3.45))  as  the  "confidence"  on 
the  initial  conditions  9(0). 

The  least  squares  algorithm,  as  can  be  seen  from  the 
simulation  results,  has  much  faster  convergence  than  the 
projection  algorithm. 

D.   PERSISTENCY  OF  EXCITATION 

In  order  to  guarantee  global  stability  of  indirect  adap- 
tive control  algorithms,  the  input  to  the  plant  has  to  be 
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persistently  exciting,  which  in  turn  implies  global  conver- 
gence of  the  plant  parameters  to  the  corresponding  actual 
values  [Ref.  3].   To  obtain  consistent  estimates  of  the  plant 
parameters,  it  is  necessary  that  the  input  signal  to  the  plant 
be  sufficiently  rich  in  frequencies,  and  excites  all  modes 
of  the  plant. 

In  a  closed  loop  set-up,  the  control  input  is  the  sum  of 
an  external  input  signal  v(t)  and  a  feedback  signal  from 
the  adaptive  controller. 

The  feedback  signal  may  in  principle  cancel  any  excitation 
contained  in  the  external  input  signal.   This  problem  and 
the  potential  for  unbounded  growth  of  the  control  and  output 
signals  have  made  the  guarantee  of  persistent  excitation  a 
difficult  problem.   Recently,  Elliot  [Ref.  8]  gave  sufficient 
conditions  which  guarantee  persistency  of  excitation.   The 
following  theorems  summarize  these  conditions. 
Theorem  4.1; 

Let  w  be  the  number  of  parameters  estimated.  In  order  to 
guarantee  global  convergence  of  the  plant  parameters  to  their 
true  values,  the  following  conditions  have  to  be  satisfied: 

1.  The  external  input  v(t)  should  consist  of  a  sum  of  2w 
sinusoids . 

2.  The  compensator  parameters  should  be  updated  each  N 
samples,  with  N  >^  lOn. 

Therefore,  in  this  thesis  report,  block  processing  is 
used  in  the  sense  that  N  data  samples  are  taken,  and  N  iter- 
ations of  the  recursive  least-squares  algorithm  are  performed 
between  control  parameters  (h,k)  updates.   To  avoid  time 
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variation  of  the  plant  dynamics  during  the  spanning  process, 
the  control  parameters  are  held  constant  during  spanning 
blocks  (of  length  N)  and  are  changed  only  between  them. 

As  we  can  see  from  Theorem  4.1,  these  two  conditions 
are  sufficient  to  guarantee  parameter  convergence  for  any 
initial  conditions  as: 


lim  0 
k+°°  K 


6,  being  the  estimate  of  9*  at  time  t,  . 

In  the  next  section,  we  will  examine  finite  time  persis- 
tency of  excitation. 

E.   FINITE  TIME  PERSISTENCY  OF  EXCITATION 

It  is  clear  that  the  persistency  of  excitation  condition 
on  the  external  input  v(t)  is  the  main  limitation  about  this 
adaptive  control  algorithm.   Recently,  in  a  report  by  Cristi 
[Ref .  15] ,  a  possible  solution  to  this  problem  has  been  given, 
by  stopping  parameter  adaptation  when  the  performances  are 
close  to  the  desired  ones.   More  precisely,  the  model  output 
error  between  desired  output  of  the  closed-loop  system  and 
controlled  plant's  output  is  the  measure  of  how  far  the  system 
is  from  the  desired  performance  (i.e.,  pole  placement)  and 
adaptation  is  stopped  whenever  the  error  falls  below  an  arbi- 
trary, preassigned  threshold.   The  configuration  of  the  system 
is  given  in  Figure  3.5. 
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e 
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l  XL 


CONTROLLER 


PLANT 


CONTROLLED 


rt(D) 

p*(D) 


f- 


v(t) 


^0 


y(t) 


Figure  3.5.   Indirect  Adaptive  Control  Scheme  with 
Added  Finite  Time  P.E. 


The  model  output  error  can  be  defined  as 


r  (D) 
u(t)   =   y(t)  -  ^^y  v(t; 


(3.48) 


The  finite  time  persistency  of  excitation  algorithm  can 
be  given  as  follows: 

1.  t  =  t+1. 

2.  Compute  \i(t)    from  Equation  (3.48). 

3.  If  | u ( t ) |  >  £/  compute  compensator  parameters  as  in 
Chapter  III. 

4.  If  |y (t) |  <  e  go  to  1. 
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A  major  difficulty  with  this  algorithm  is  to  be  able  to 
guarantee  global  stability  of  the  whole  system.   Therefore 
we  have  to  prove  that  the  signals  in  the  loop  are  bounded, 
provided  the  external  input  v(t)  is  bounded.   Global  stability 
of  the  closed-loop  system  is  proved  below. 

When  time  goes  to  infinity  and  persistency  of  excitation 

is  present,  parameters  of  the  estimates  of  the  plant  poly- 

nomials  r.  and  P   converge  to  the  actual  values  r  and  p. 

Therefore,  the  difference  between  measured  output  and  desired 

output  tends  to  zero,  and  also  U(t)  tends  to  zero.   Thus,  by 

definition  of  limit,  there  exists  an  instant  t^  such  that 

o 


u(t) 


for  all  t  >  t   and  for  any  given  e  >  0.   Rearranging  Equation 


(3.48)  yields 


y(t)   = 


_   r(D) 


p*(D) 


v(t)  +  y(t) 


(3.49) 


which  can  be  expressed  on  a  block  diagram  form  as  in  Figure  3.6 


Figure  3.6. 


Block  Diagram  Representation  of  Equation 
(3.49) 
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It  can  be  interpreted  as  a  single  input  single  output 
linear  stable  process  with  bounded  input  v(t)  and  bounded 
disturbance  y (t) .   From  the  fact  that  stable  linear  systems 
with  a  bounded  input  produce  a  bounded  output,  the  global 
stability  results  follow  easily.   So  that  with  finite  time 
persistency  excitation,  the  closed-loop  system  is  eventually 
equivalent  to  a  system  with  poles  as  desired,  bounded  zeroes 
with  a  bounded  output,  disturbance  y  (t)  . 

We  can  generate  the  external  input  v(t)  as 


v(t)   =   ve(t)  +  vc(t)  (3.50) 


with  v  (t)  the  external  desired  command  and  v  (t)  an  added 
c  e 

persistency  excitation  signal.   When  the  adaptation  is  con- 
tinuing, the  external  input  is  composed  of  v  (t)  and  v  (t) , 
but  after  stopping  the  adaptation  persistency  of  excitation 
is  not  needed,  it  will  be  identical  to  the  longer,  so  that 
v  (t)  can  decay  to  zero. 

In  this  way,  the  adaptive  algorithm  is  activated  only 
when  the  performance  index  y(t)  is  larger  than  a  minimum 
threshold.   In  the  next  chapter,  simulation  studies  will  show 
their  efficiency  via  some  examples. 
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IV .   SIMULATION  STUDIES 

This  research  report  includes  three  computer  programs. 
The  program  named  CONDIS ,  given  in  Appendix  B,  is  used  to 
investigate  the  behavior  of  the  zero  in  sampled  data  systems. 
This  program  computes  the  parameters  of  the  discrete  time 
transfer  function  from  the  parameters  of  its  continuous  time, 
and  the  sampling  rate.   This  program  has  been  used  to  inves- 
tigate the  perturbations  of  the  zeroes  for  different  values 
of  the  sampling  interval. 

The  other  two  programs  given  in  Appendix  C  and  D  simulate 
the  indirect  adaptive  control  using  recursive  least  squares 
and  projection  algorithms.   Entering  order  and  numerator, 
denominator  parameters  of  the  plant,  observer  polynomial  and 
desired  closed-loop  characteristic  polynomial,  the  program 
simulates  the  entire  system  in  an  interactive  fashion.   Also 
it  is  capable  of  graphical  and  tabulation  results. 

The  adaptive  controller  presented  above  has  been  simulated 
for  several  different  plant  dynamics. 
Example  4.1; 

Let  the  discrete  time  transfer  function  of  the  plant  be: 


H(z)   =      Z    +    2  2  +  2 


z2-2z+0.75      (2  -1-5)  U  -0.5) 


The  plant  has  one  unstable  zero,  z  =  -2,  and  two  poles 
p,  =  1.5  and  p„  =  0.5. 

50 


The  polynomials  q(z)  and  p*(z)  are  chosen  to  be  stable 
polynomials  with  degree  n. 
In  particular,  let 


q(z)   =   z2  -  0.3z  -  0.28   =   (z  -0.7)  (z  +0.4) 


and 


p*(z)   =   z2  -  l.lz  +  0.3   =   (z  -0.5)  (z  -0.6) 


both  having  stable  roots,  i.e.,  inside  the  unit  disc  in  the 
z-plane.   Before  going  into  simulation  study,  the  problem  is 
solved  analytically  by  comparing  their  responses.   By  writing 
the  plant  dynamics  in  shift-operator  form, 


-1 


"I  _>.->    "2 


H(q  *)   =        _-,        _2 
1  -2q    +0.75q 

we  can  derive  the  difference  equation  of  the  given  system  as 


y(t)   =   2y(t-l)  -  0.75y(t-2)  +  u(t-l)  +  2u(t-2) 


y(t)   =   $(t-l)T6* 


where 
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and 


(t-1 


-y(t-l) 

-y(t-2) 

u(t-l) 

u(t-2) 


* 


-2 
0.75 
1 
2 


9*  being  the  vector  of  the  actual  parameters  of  the  plant 
Generally,  it  can  be  partitioned  as: 


6*   = 


p  and  r  being  vectors  of  parameters  of  the  denominator  and 
numerator  of  the  plant,  respectively. 

From  Equation  (3.23),  considering  k(z)  and  h(z)  as 
controller  polynomials,  we  can  write 

k(z)p(z)  +  h(z)r(z)   =   q(z)[p(z)  -p*(z)] 

In  this  example  we  have  chosen: 
p(z)   =   z2  -  2z  +  0.75 


52 


r(z)   =   z  +  2 

q(z)   =   z2  -  0.3z  -  0.28 

p*(z)   =   z2  -  l.lz  +  0.3 

By  the  constraints  on  the  polynomials  for  solvability  of 
the  Diophantine  equation,  we  choose  k  and  h  to  be  first  order 
polynomials  as: 


k(z)   =   k,z  +  k 


h(z)   =   h,z  +  h 


Substituting  these  polynomials  into  the  Diophantine 
equation  in  matrix  form  yields: 


0.75   -  2    1    0 


-0.126 


10    0 


0.117 


0     0.75  -2    1 


0.72 


2    10 


-0.9 
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Since  the  plan  does  not  have  any  pole-zero  cancellation, 

the  determinant  of  the  matrix  S   is  not  equal  to  zero.   Hence, 

m  A 

it  is  solvable.   Solving  the  above  matrix  equation,  we  compute 
the  control  parameters  k,  =  -.899,  k   =  -0.688,  h,  =  -0.3900 
and  h   =  0.195.   The  corresponding  polynomials  of  the  con- 
troller are: 


k(z)       =•     -.899z    -    0.688 
h(z)       =      0.39z    +    0.195 

This  corresponds  to  the  optimal  choice  of  compensator 
parameters  for  the  desired  pole  placement.  Therefore  the 
compensator  given  by  the  difference  equation 


=   MD)u(t)  +  h<D)y(t)  + 
q  (D) 


yields  closed-loop  transfer  function 


H(z)   = 


Y(z) 
V(z) 

q  (z)  r  (z) 


p(z)  [q(z)  -k(z)]  -r(z)h(z) 


which  becomes: 


H(z)   = 


z  +2 


2 
z   -  1  .lz  +  0 . 3 
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The  step  response  of  this  system  is  given  in  Figure  4.1. 
The  next  approach  to  the  problem  is  simulation  of  the  system 
using  a  computer  program.   It  is  observed  that  from  the  simu- 
lations, after  the  lOn  iterations,  the  estimated  parameter 
will  be  closer  to  6*  and  6  will  tend  to  zero,  where 


*  _ 


This  means  that  p  and  r  are  converging  to  the  actual  param- 
eters,  and  h  and  k  converge  to  the  actual  values  for  the 
desired  pole-placement.   The  desired  output  and  controlled 
system's  output  are  given  in  Figure  4.2.   Also,  parameter 
error  and  output  prediction  error  are  given  in  Figures  4.3 
and  4.4,  respectively.   After  convergence  of  the  prediction 
error  to  zero,  the  persistency  of  excitation  added  to  the 
external  command  is  turned  off.   In  this  example,  we  have 
chosen  the  step  as  external  command.   External  input  and  plant 
input  are  in  Figures  4.5  and  4.6,  respectively. 

If  we  compare  Figures  4.1  and  4.2,  analytical  and  simu- 
lated system's  results  are  almost  the  same. 

In  the  next  example,  it  is  considered  that  one  of  the 
plant  parameter  is  changed  after  some  period  of  time.   It  is 
investigated  how  the  system  behavior  will  change  in  this 
condition. 

Example  4.2: 

Let  the  plant  be  the  same  as  the  one  in  the  previous 
example.   After  two  blocks  length,  the  zero  of  the  plant  is 
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perturbed  from  z  =  -2  to  z  =  -3.   To  compute  the  parameters 
associated  to  the  perturbed  plant  with  pulse  transfer  function 


H(z)   = 


z  +  3 


z   -  2z  +  0.75 


write  its  difference  equation  as 


y(t)   =   2y(t-l)  -  0.75y(t-l)  +  u(t-l)  +  3u(t-2) 


and  immediately 


* 


-2.0 
0.75 
1.0 
3.0 


By  using  the  same  observer  and  desired  polynomials  g(z), 

p*(z)  and  following  the  same  steps  as  in  Example  4.1,  the 

controller  parameters  become  k,  =  -0.894,  k   =  -0.778, 
r  1  o 

hn  =  -0.305  and  h   =  0.152. 
1  o 

Then,  the  polynomials  k(z)  and  h(z)  are: 


k(z) 


!94z  -  0.778 


h(z)   =   -0.305z  +  0.152 


The  closed  loop  transfer  function  becomes 
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H(z)   =    2   Z  +  3 

z   -  l.lz  +  0.3 


The  step  response  of  this  transfer  function  is  given  in 
Figure  4.7.   Also,  the  other  graphical  results  are  in  Figures 
4.8-4.12. 

In  the  next  example,  we  investigate  how  the  disturbance 
at  the  output  will  affect  the  behavior  of  the  controlled 
system. 

Example  4.3; 

Consider  a  plant  with  pulse  transfer  function  as: 


H(»   = 


z2  -  2z  +  0.75 


and  assume  an  output  disturbance  exists.   Let  the  disturbance 
be  sinusoidal  with  frequency  5tt/2  rad/sec.   In  this  case,  it 
is  observed  that  the  estimated  parameters  of  the  plant  don't 
converge  to  the  actual  parameters.   Corresponding  plots  are 
given  in  Figure  4.13  through  4.17.   However,  if  the  disturbance 
is  small,  the  parameters  converge  close  to  the  actual  values 
and  we  can  still  obtain  satisfactory  performances. 

Comparison  of  Different  Estimation  Techniques:   RLS  and  P. A. 

In  the  above  examples  we  used  a  recursive  least-squares 
algorithm  for  estimating  the  plant  parameters.   In  this  sec- 
tion, we  compare  the  behavior  of  the  indirect  adaptive  control, 
when  the  projection  algorithm  is  used. 

The  projection  algorithm  has  the  advantage  of  requiring 
less  computations.   Approximately,  the  number  of  computations 
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2 
grows  as  n   in  the  recursive  least-squares  algorithm.   But 

in  the  projection  algorithm,  it  grows  as  n.   As  noted  before, 

the  projection  algorithm  can  be  described  as 


(t)   =   0(t-l)  +  a(Mt  ^ [y(t)  -$(t-l)T9  (t-1)] 

c  +$(t-l)  $(t-l) 


with  c  >  0  and  0  <  a  <  2.  The  value  of  the  parameter  a  is 
crucial.  The  behavior  os  the  algorithm  greatly  depends  on 
its  value. 

In  the  next  example,  the  projection  algorithm  is  used 
for  estimation  procedure  on  the  previous  process. 
Example  4.4: 

Using  the  same  pulse  transfer  function  and  controller 
polynomials  as  in  Example  4.3  and  choosing  the  constants 
a  =  1  and  c  =  1  in  the  above  equation,  results  can  be  observed 
from  Figure  4.18  through  4.22. 

Actually,  using  this  algorithm,  several  systems  have 
been  simulated.   This  result  is  the  most  reasonable  one. 

Considering  the  above  results,  indirect  adaptive  control 
is  a  very  effective  control  technique  when  the  parameters  of 
the  plant  are  unknown  and  large  and  unpredictable  variations 
occur  in  the  plant  dynamics.   Also,  it  is  applicable  to 
nonminimum-phase  systems.   This  can  be  considered  as  a  big 
advantage  of  the  method. 
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APPENDIX  A 
DESCRIPTION  OF  COMPUTER  PROGRAMS 

Computer  programs  are  written  using  WATFIV  and  FORTRAN 
programming  languages.   All  of  them  are  prepared  as  an 
interactive  program.   The  following  explanation  is  about  how 
one  can  use  these  programs . 

The  program  named  CONDIS  given  in  Appendix  B  finds  the 
pulse  transfer  function  from  the  continuous  time  transfer 
function.   The  program  asks  for  the  sampling  interval,  orders 
and  coefficients  of  the  continuous  time  transfer  function 
polynomials.   The  continuous  time  system  is  described  in  state 
variable  form  as 

x   =   Ax  +  Bu 

and  the  corresponding  discrete  time  system  as 

x(k+l)   =   $x(k)  +  Tu(k) 


This  program  gives  the  $  and  V   matrices  corresponding  to  A 
and  B. 

There  is  no  limitation  about  the  order  of  the  system  i 
all  programs.  For  higher  order  systems,  the  dimensions  of 
the  declared  matrices  should  be  increased. 

The  program  named  RCIOP  given  in  Appendix  D  simulates 
an  indirect  adaptive  control  system  using  the  projection 


n 


algorithm.   The  program  asks  some  questions  to  the  user  in 
the  following  order: 

1.  The  order  of  a  plant  (N)  ; 

2.  The  constants  a  and  c  that  are  given  in  Equation  (3.38); 

3.  Coefficients  of  the  numerator  polynomial  of  the  plant 

(in  ascending  order  of  z) ; 

4.  Coefficients  of  the  denominator  polynial  of  the 
plant  (in  ascending  order  of  z)  ; 

5.  Coefficients  of  the  controller  polynomials  g(z)  and 
p*(z) . 

The  external  input  is  defined  inside  the  program.   It 
can  be  changed  when  it  is  necessary. 

The  last  program  RCIOR  given  in  Appendix  C  is  the  most 
important  one.   Least-squares  algorithms  are  used  for  parameter 
estimation.   This  program  asks  the  same  questions  about  system, 
except  question  #2. 

Computer  programs  RCIOR  and  RCIOP  give  the  graphical 
and  tabulated  results.   Graphical  outputs  are  obtained  using 
DISSPLA. 
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APPENDIX  B 
COMPUTER  PROGRAM  CONDIS 


$  JOB    WATF I V   ONUK , XRE F , EXT 

C        DEFINITION  OF  VARIABLES: 

INTEGER  N,M,R,  I  ,K,L,  Ll/L2/  II,  J,  I3,M1,M2,N3 

integer  M3 , L3 , 16 

REAL   AA( 2 , 2 ) , BB( 2 , 1 ) , ID( 2 , 2 ) , PS( 10 , 10 ) 

REAL   FI    i6/I6)/GT(10/10)/61( iO/10)/C2( 10,10) 

real    c3( 10, 10) ,081 10,10),d4    10,10),d9( 10,10 

REAL   C4£10,  10)  ,C5(  10,  10)  ,C6(  10,  10)  ,C7(  10,  10) 

REAL   C10( 10 

REAL 

real    cl4(  I0,l0),cl9(  I0,10),c20(  10 ;  10) )  e4(_  10  i  10) 

REAL   D10( 10,10    ,D11( 10,10    ,C12( 10,10),C13( 10,10) 

REAL   C15(10,10),C16(  10,10    ,C17(  10,  10)  ,  018^10,  10) 


,L   C4C10,10),C5[10,10),C6( 10,10), C7(  10.10) 

,l  cio(  16, 16  ,cli(io. 16) ,D2( io, 16) ,D3(io. io) 

,L  D5(  10,.l0),D6(10,l6),D7(  10,l6),D8(  10,  l6) 
.1    cl4( l6,l6j,cl9( i0,l6),c20(  10,l0),e4(l0,10 
,L  D10( 10,10    ,D11( 10,10    ,C12( 10,10),C13( 10 

»d*  cisfio^ioLciei  10,10  ,C17(  10,10), cistio 

REAL  C110(10, 10) ,W1(1, 10)/l0*0. /, ID1( 10, 10) 

REAL  CC( 10),TT( 16) ,DD[ 10  , RR( 2 , 2) , PP(2, 2 ) 

REAL  Dl( 10/l0),El(i0,10),E2(  10, l6) , E3( 10, 10 ) 

REAL  E5  10,10i,E6( 10, 10  j , C9[ 10, 10), ss( 2 , 2 ) 

REAL  FF(2,1),GG( 1 , 1 ) , ALPHA( 10) , CO( 10 ) , ALP, TAR 

REAL  TC . TRG, TR1 , TR2 , TR3 , TR4, TR5 . TR6 , TR7 , TR8 

real  TR9, TRlO, TR11, TK1, c21(l0,  l6 ) 

PRINT/ THIS  PROGRAM  FINDS  THE  Z-DOMAIN  PULSE  TRANSFER' 

PRINT, 'FUNCTION  WHEN  GIVEN  THE  S-DOMAIN  TRANSFER  FUNC. 

PRINT, '  T 

PRINT, 'THE  S-DOMAIN  TRANSFER  FUNCTION  SHOULD  BE 
IN  THIS  FORM' 
C 
C       H(S)=NUM(S)/DEN(S) 

C  NUM(  S)=T(M)S**N-1+ +T(2)S+T(  1) 

C  DEN( S)=S**N+D(N)S**N-1+ +D(  2)S+D(  1) 

PRINT, '  ' 

PRINT, 'ENTER  THE  NUMERATOR  POLYNOMIAL  DEGREE' 

READ, Ml 

PRINT, T  ' 

PRINT, 'ENTER  THE  DENOMINATOR  POLYNOMIAL  DEGREE' 

READ,N 

M2=Ml+l 

DO  2  I1=1,M2 

PRINT,  fT(  '  , II,  '  )=?' 
READ , TT( I 1 ) 
2       CONTINUE 

N3=N-1 

M3=Ml+2 

IF  (Ml.  LT.N3)  THEN 
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DO  82  I1=M3,N 
TT( I1)=0. 0 
82         CONTINUE 
END  IF 


DO  12  H=1,N, 

PRINT, TD(  ' , II, ' )=?' 
READ , DD( I 1 ) 
12      CONTINUE 

DO  22  I1=1,N3 
J=l 

AA( II, J)=0. 0 
22      CONTINUE 

DO  32  J=2,N 

DO  102  I1=1,N3 
13=11+1 
IF  (J. EQ. 13)  THEN 

AA(  II, J)=l. 0 
ELSE  * 

AA( II, J)=0. 0 
END  IF 
102        CONTINUE 
32      CONTINUE 

DO  72  J=1,N 

AA(N, J)=-DD( J) 
72      CONTINUE 

DO  52  I1=1,N3 
J=l 

BBf II, J)=0. 0 
52      CONTINUE 

BB(N£J)=1.  0 
DO  62  I1=1,N 

CC(  I1)=TT(  II) 
62      CONTINUE 
R=l 

M=N  ,   , 

CALL  FAMILY(N,M, 'AA' ) 
CALL  RESULT(N,M, I.K.AA) 
CALL  FAMILY(N,R, TBB* ) 
CALL  RESULT  N,R, II, L3,BB) 
CALL  FAMILY( R,N, TCCT) 


:i* 


CALL  RESULT(R,N, I,K,CC) 
DO  14  J=1,N 

DO  15  1=1 -N 

IF  ( I. EQ. J)  THEN 

ID( I, J)=l. 0 
ELSE 

ID(  I, J)=0.  0 
END  IF 
15         CONTINUE 
14      CONTINUE 
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CALL  FAMILY(N,M, ' ID' ) 
CALL  RESULT) N,M, I, J, ID) 
68        M=l 

DO  285  1  =  1  ,N 

DO  286  J=1,N 

RR( I, J)=ID( I, J) 
286  CONTINUE 

285       CONTINUE 
TR1=1. 0 

CALL  ALI(N.N,  I,K,AA,W1) 

GAMMA=AMAX1[W1(1, 1)  .Wl(  1.2)  .Wl(  1,3)  .Wl(  1,  4) 
*     wl(l,  5}a{J1(-L'6)  Mil,  7)  ,fil(i,  8)  ,ftl(i,  9)  ,fil(l, 10)) 

PRINT.TC 

THl=TRl/2 

PRINT, TH1 

TH2=TRl/4 

PRINT,  TH2- 

TH3=TRl/8 

PRINT.TH3 

TH4=TR1/16 

PRINT, TH4 

IF  (TR1.  LE.TC)   THEN 

TR2=TRl/2 

TR3=TRl/3 

TR4=TRl/4 

TR5=TRl/5 

TR6=TRl/6 

TR7=TRl/7 

TR8=TRl/8 

TR9=TRl/9 

TR10=TR1/10 

TR11=TR1/11 

CALL  SCALAR(N,N, I,K,TR1, ID, CI) 

CALL  SCALAR) N,N, I ,K,TR2, AA,D2 ) 

CALL  CALCUL  N,N,N, I , K, K, C2 , CI ,D2 ) 

CALL  SCALAR) N, N, I, K,TR3,AA,D3) 

CALL  CALCUL(N,N,N, I , K, K, C3 , C2 ,D3 ) 

CALL  SCALAR) N, N, I, K,TR4,AA,D4) 

CALL  CALCUL(N,N,N, I , K, K, C4, C3 ,D4) 

CALL  SCALAR  N,N, I, K,TR5,AA,D5) 

CALL  CALCUL(N,N,N, I , K, K, C5 , C4, D5 ) 

CALL  SCALAR) N, N, I, K,TR6,AA,D6) 

CALL  CALCUL(N,N,N, I , K, K, C6, C5 , D6 ) 

CALL  SCALAR) N, N, I, K,TR7,AA,D7) 

CALL  CALCUL(N,N,N, I , K, K, C7, C6 . D7 ) 

CALL  SCALAR  N,N, I, K,TR8,AA,D8) 

CALL  CALCUL)n,N,N, I , K, K, C8, C7 , D8 ) 

CALL  SCALAR) N, N, I, K,TR9,AA,D9) 

CALL  CALCUL)N,N,N, I , K, K, C9, C8, D9 ) 
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CALL  SCALAR(N,N, I,K,TR10,AA,D10) 
CALL  CALCUL(N,N,N,  I ,  K,  K,  CIO,  C9.D10  ) 
CALL  SCALAR  N,N, I, K,TRll,AA,Dll) 
CALL  CALCUL  N.N.N, I , K. K, Cll , CIO, Dll ) 
CALL  SUM(N,N, I , K, CI . C2 , C12) 
CALL  SUM(N,N, I , K, C12 , C3 , C13 
CALL  SUM  N,N, I , K, C13 , C4, C14 
CALL  SUM  N,N, I , K, C14, C5 , C15 
CALL  SUM(N,N, I , K, C15 , C6 , C16 
CALL  SUM(N,N, I , K, C16, C7, C17 
CALL  SUM  N,N, I , K, C17 , C8, C18 
CALL  SUM  N,N, I , K, C18, C9 , C19 . 
CALL  SUM  N,N, I , K, C19 . CIO, CllO) 
CALL  SUM(N,N, I, K, CllO, Cll. PS) 
CALL  CALCUL( N, N. N, I , K . K, C21 , AA, PS ) 
CALL  SUM(N,N, I,K. ID.C21,FI) 
CALL  FAMILY(N,N,  'PSM 
CALL  RESULT  N,N, I.K.PS) 
CALL  FAMILY) N, N, TFIf) 
CALL  RESULT  N,N, I, K,FI) 
CALL  CALCUL  N,R,N, I , K, K, GT, PS, BB) 
CALL  FAMILY(N,R,  GT' ) 
CALL  RESULT(N,R, I,L,GT) 
ELSE 
TRG=TR1 
L1=0. 0 
TK2=2  0 
677       TRG=TR1/(2**L1) 

IF  (TC. LE. TRG)  THEN 

L1=L1+1 

GO  TO  677 

END  IF 

PRINT. LI 

TR=TR1/(2**L1) 

CALL  SCALAR(N,N, I,K,TR, ID, CI) 

TR2=(TR**2)/2 

CALL  SCALAR(N,N, I , K, TR2 , AA, C2 ) 

TR3=(TR**3)/6 

CALL  CALCUL(N,N,N, I , K, K,D1, AA, AA) 

CALL  SCALAR(N,N, I , K, TR3 ,D1, C3 ) 

TR4=(TR**4)/24 

CALL  CALCUL(N,N,N, I , K, K,D2 , AA, Dl ) 

CALL  SCALAR  N.N, I, K,TR4,D2,C4) 

TR5=(TR**5)/120 

CALL  CALCUL(N,N,N, I , K, K, D3 , AA, D2 ) 

CALL  SCALAR  N.N, I, K,TR5,D3,C5) 

TR6=(TR**6)/720 

CALL  CALCUL(N,N,N, I , K, K, D4, AA, D3 ) 

CALL  SCALAR  N.N.I, K,TR6,D4,C6) 

TR7=(TR**7)/5040 
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CALL  CALCUL(N/N/N/ I , K, K,D5 , AA, D4) 

CALL  SCALAR) N. N. I, K,TR7,D5,C7) 

TR8=(TR**8)/40320 

CALL  CALCUL(N,N,N,  I ,  K,  K,  D6,  AA.D5  ) 

CALL  SCALAR) N,N, I, K, TR8,D6,C8 ) 

TR9=( TR**9 )/362880 

CALL  CALCUL(N,N,N, I , K, K,D7, AA, D6) 

CALL  SCALARtN,N, I.K,TR9,D7,C9) 

TR10=£TR**10)/3628800 

CALL  CALCUL(N,N,N, I , K, K.D8. AA. D7 ) 

CALL  SCALAR) N,N. I, K, TR10,D8, ClO) 

TRll=(TR**ll)/399l6800 

CALL  CALCUL(N,N,N, I , K, K.D9 . AA.D8) 

CALL  SCALAR  N,N, I, K, TRll, D9, Cll ) 

CALL  SUM(N/N,I,K,C1.C2ZC12] 

CALL  SUM)N,N, I , K, C12 , C3 , C13 

CALL  SUM)N,N, I , K, C13 , C4, C14 

CALL  SUM  N,N, I , K, C14, C5 , C15 

CALL  SUM(N,N, I , K, C15 , C6, C16 

CALL  SUM)n,N, I , K, C16, C7, C17 

CALL  SUM  N,N, I , K, C17, C8, C18 

CALL  SUM  N,N, I , K, C18, C9 . C19 , 

CALL  SUM(N,N, I , K, C19 , CIO, C20) 

CALL  SUM)N,N,  I , K. C20, Cll, PS ) 

CALL  FAMILY(N,N,  fPlM 

CALL  RESULT) N, N, I, K, PS) 

DO  255  16=1, Ll 

CALL  SCALAR(N,N,I,K.TK2, ID, El) 

CALL  FAMILY) N,N,  El* ) 

CALL  RESULT) N,N, I,K,E1) 

CALL  CALCUL)N,N,N, I , K, K, E2 , AA, PS) 

CALL  SWlTHjNjIjK.tl.E^,^) 

CALL  FAMILY(N,N, fE3* ) 

CALL  RESULT) N,N, I, K,E3) 

K,K,E5,E3,PS) 


:,£5) 
CALL  COPY(N,N, I,K,PS,E5) 
CALL  FAMILY) N, N, ' PST j 
CALL  RESULT  N,N, I,K,PS) 
255       CONTINUE 


CALL  FAMILY(N,N, 'PS' ) 

~",PS) 

",K,GT,PS,BB) 


CALL  RESULT) N,N, I ,K, 


T) 

CALL  CALCUL)  N'U'  fH'l ',  K,  K,  E6 ,  AA,  PS  ) 
CALL  SUM(N,N, I,K. ID,E6,FI) 
CALL  FAMILY(N,N, 'fI* ) 
CALL  RESULT  N,N, I, K,FI) 
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END  IF 

CALL  CALCUL(N/M/N/ I/K/K/FF/RR/GT) 

CALL  CALCUL(M,M,N, I , K, K, GG, CC, FF) 

CO(N)=GG(M,M) 

PRINT,   > 

PRINT, '  ' 

PRINT, "NUMERATOR  COEFF. S  OF  THE  PULSE  TRANSFER 

FUNCTION1 
PRINT, '  T 

PRINT , ' NUMCOEF ' , CO( N ) 
DO  147  L=1,N 
EXECUTE  AR 
ALP=-TAR/L 

CALL  CALCUL(N,N,N, I, J, J,PP,FI,RR) 
CALL  SCALAR(N,N, I, J, ALP, ID,SS) 
CALL  SUM(N,N, I,J,PP,SS,RR) 
CALL  CALCUL(N,M,N, I , J, J, FF, RR, GT) 
CALL  CALCUL(M,M,N, I , J, J, GG, CC, FF) 
CO(  L)=GG(M,M) 
ALPHA(  L ) =ALP 

147  CONTINUE 

DO  150  L=1,N3 

PRINT, TNUMCOEF* ,CO(L) 

150  CONTINUE 
PRINT, r 
PRINT, '  ' 

PRINT, 'DENOMINATOR  COEFF. S  OF  THE  PULSE  TRANSFER 

FUNCTIONT 
PRINT. '  T 
DO  151  L=1,N 

PRINT, TDENUMCO' ,ALPHA(  L) 

151  CONTINUE 
STOP 

REMOTE  BLOCK  AR 

TAR=0. 0 

CALL  CALCUL(N,N,N, I , K, K, PP, FI , RR) 

DO  148  1=1, N 

DO  149  J=1,N 

IF(I.EQ.  J)  THEN 

TAR=TAR+PP(  I, J) 
END  IF 
149  CONTINUE 

148  CONTINUE 
END  BLOCK 
END 

C   ***THis  SUBROUTINE  READS  THE  MATRICES  FROM  THE  DATA 
FILE. *** 
SUBROUTINE  ENTER( J, D, E, F, G) 
INTEGER   J  ,  D  ,  E  ,  F 
REAL  G( J,D) 
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DO  60  E=1,J 

READ  50.(G(E,F),F=1,D) 
50        FORMAT( 10F8. 4) 
60      CONTINUE 
RETURN 
END 
C   ***THIS  SUBROUTINE  PRINTS  THE  ARRAYS  AS  AN  MATRIX 
FORM. *** 
SUBROUTINE  COPY( H, 0, P, S, Gl , Kl ) 
INTEGER  H,0,P. S 
REAL  Gl(H,OJ  ,K1(H,0) 
DO  131  P=i,H 

DO  132  £=1,0 

G1(P,S)=K1(P,S) 
132        CONTINUE 
131     CONTINUE 
RETURN 
END 

SUBROUTINE  RESULT( H, 0, P, S, T) 
INTEGER  H  ,    0  ,  P  ,  S 
REAL  T  (  H  .  0  ) 
DO  80   P  =  1  ,  H 

PRINT  70,(T(P/S).S=1.0J 
70        FORMAT(  T6> , T2 , lOX, 10( 3X, F12. 6)) 
80      CONTINUE 
RETURN 
END 
C   ***THIS  SUBROUTINE  CALCULATES  THE  MULTIPLICATION  OF 
TWO  MATRICES. *** 
SUBROUTINE  CALCUL( U. V, Y, Z,X, ZX, ZT, W, Q) 
INTEGER  U,V,Y,Z,X,ZX 
REAL  ZT(U/V)/W(U/Y)/Q(Y/V) 
DO  93  Z=1,U 
DO  92  X=1,V 

ZT(Z.X)  =  0.  0 
DO  91  ZX=1,Y 

ZT( Z,X)=ZT( Z,X)+W( Z,ZX)*Q( ZX,X) 

91  CONTINUE 

92  CONTINUE 

93  CONTINUE 
RETURN 
END 

C  ***THIS  SUBROUTINE  WRITES  THE  MATRIX  NAME  AND  ROW, 
COLUMN  NUMBERS. ***/ 

SUBROUT I NE  FAM I LY(  ROW , COLUMN , NAME ) 
INTEGER  ROW, COLUMN 
CHARACTER* 2  NAME 
PRINT  94. 'MATRIX*  NAME 

94  FORMAT(Tl. '0T,T2, lOX, T12 . A6, T19, IX, T20, A2 ) 
PRINT  95/ ROW  NUMBER.  T  ,  ROW 
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95  FORMATS Tl. ' 0 '  T2 , 10X, T12 , All , T23 , 8X, T31, 12) 
PRINT  96. f COLUMN  NUMBER.  f . COLUMN 

96  FORMAT(Tl,T0r,T2,10X,T12,A14,T26,5X,T31, 12) 
RETURN 

END 
C   ***THIS  SUBROUTINE  CALCULATES  THE  SUMMATION  OF  THE 
C      absolute  values  OF  THE  SAME  COLUMN  ELEMENTS. *** 
SUBROUTINE  ALI( II , 12 ,H1 ,H2 , Gl, W2 ) 
INTEGER  II, 12, HI, H2 
REAL  W2(l, 12), Gil  11,12) 
DO  41  H2=l, li 

DO  42  61=1,11 
W2(  1,H2)=W2(  1,H2)+ABS(G1(H1,H2) ) 

42  CONTINUE 
41      CONTINUE 

RETURN 
END 
C   ***THIS  SUBROUTINE  MULTIPLIES  THE  MATRIX  BY  SCALAR 
NUMBER. *** 
SUBROUTINE  SCALAR( Kl , K2 , Gil , G2 , PI , G3 , G4 ) 
INTEGER  K1,K2,G11.G2 
REAL  G3(K1,K2_)/LG4(K1/K2)/P1 
DO  43  G2=1,K2 

DO  44  G11=1,K1 
G4(G11/G2)=G3(G11/G2)*P1 

44  CONTINUE 

43  CONTINUE 
RETURN 
END 

C   ***THIS  SUBROUTINE  CALCULATES  THE  SUM  OF  THE  TWO 
MATRICES. *** 
SUBROUTINE  SUM(K3.K4, 15, K5,X1,X2 ,X3 ) 
INTEGER  K3,K4/I5.k5 

REAL  Xl( K3 , K4) , X2( K3 , K4) , X3(  K3 , K4) 
DO  45  K5=1,K4 

DO  46  15=1. K3 

X3( I5,K5)=X1( I5,K5)+X2(  I5,K5) 
46      CONTINUE 

45  CONTINUE 
RETURN 
END 

$ ENTRY 
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APPENDIX  C 
COMPUTER  PROGRAM  RCIOR 


rj  ****************************************** 

C  *    THESIS  PROGRAM  * 

C  *    INDIRECT  ADAPTIVE  CONTROL  * 

C  *    PARAMETER  ESTIMATION  * 

C  *    USING  REC.  LEAST  SQUARE  ALGORITHM     * 

C  *  * 

n  ****************************************** 

C 

rj    ****************************************** 

C   *     DEFINITION  OF  VARIABLES:  * 

rj    ****************************************** 

INTEGER  N1,N,L,K,L1,  I ,  M,  L2  ,  J,N3  ,N2  ,  Jl ,  J2  ,  J3  ,M1 
INTEGER  N4, N5 , N6 . COLUM, J5 , M5 , J7 , J8 , J9 , L5 , M6 , L7 
INTEGER  DEG1 , DEG2 , DEG3 . K9 , M7 , L8 , M2 , M3 , L4 , M9 , M8 


V2( 1.1 , V3(4. 1),V4( 1. 1) , V5  4. lj 
Y(200)  ,  6(200).  Fl(  4,1).  FT(  1,4),P2(4,4) 

9( 4 , 4 i , vp( 266 ) , se  806 , 4 \ 


REAL  Vl(4,l 
REAL  U2{ 4,1 
real  v6(  1, 1 

REAL  TSM.tsp;TSFrTSNlPllv(266l,'VM(266")7MU(56) 
REAL  A(  20,20) . H( 20, 20) , B( 20) , C( 20, 20) ,X(20) 
REAL  D(10)  .T(20)  .  F(  20)  ,G(20J  ,  Sl(  10),S2(  10)  ,03(10) 
REAL  U3(  4, 1)  ,C5(  10)  ,C6(  10) ,MK(  50),U4(4, 1),C7( 100) 
REAL  SUM2,TEMP,SIGMA.ASIGMA,DY,UP( 200), Q2( 50) 
REAL  Ql(50)/Q3(50),C8(_100),R(26),kD(  l66j,VC(200) 
REAL  Bl(20)/Hlf 100)lH2flOO),H3CiOO)/H4(  100) 
REAL  Ul[4.i),P3(4,l),AB(4,4),C4( 10  ,h5( 100) 
PI=3. 141592654 
DO  7777  1=1,2 

Hl( I)=0. 6 

H2  I  =0. 0 

H3  I  =0- 0 
H4  I)=0. 0 
H5( I 1=0. 0 
7777   CONTINUE 

DO  20  I  =1,200 
DO  8   J  =1,4 
SE(  I, J)=0. 0 
8       CONTINUE 
20    CONTINUE 
M=l 

Q  ******************************************** 

C   *  THIS  PART  OF  THE  PROGRAM  ASKS  THE  SYSTEM  * 
C   *  ORDER  ,  PLANT  AND  MODEL  NUMERATOR, DENOM I. * 
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C   *  POLYNOMIAL  COEFFICIENTS.  * 

n      ******************************************** 

WRITE(6.701) 
701    FORMAT( flT, 'ENTER  THE  ORDER  OF  THE  SYSTEM') 

READ( 5 ,  * )  N 

N1=N+1 

N2=2*N 

N3=N+2 

N6=2*N+1 
C      M7=10*N*2 

L8=10*N2 

CALL  FRTCMS(_  CLRSCRN  ') 

WRITEfS^OO?) 

2002  FORMAT(  '  T,TENTER  PLANT  NUMERATOR  POLYNOMIAL 

COEFFICIENTS1) 
WRITE( 6.2003) 

2003  FORMAT [ f  T.TIN  ASCENDING  ORDER  OF  Z ' ) 
DO  2001  1  =  1, N 

J7=I-1 
WRITE(6.2004)  J7 

2004  FORMAT (  i    T,rCOEFF.  OF  Z**(  ' , 12 ,  '  )  =  '  ) 
READ(5>)  C5(  I) 

2001   CONTINUE 

CALL  FRTCMS(  ' CLRSCRN  '  ) 
WRITE(6.2005) 

2005  FORMAT( i    T, T ENTER  PLANT  DENOMINATOR  POLYNOMIAL 
COEFFICIENTS1) 


WRITE(6,2006) 
FORMATC    ,  IN 


2006  FORMATS  T.rIN  ASCENDING  ORDER  OF  Z') 
WRITE(  6.2007) 

2007  FORMATS '  T. 'HIGHEST  COEFF.  SHOULD  BE  1.  ') 
DO  2008  I  =  1,N 

J8=I-1 

WRITE(6.2009)  J8 
2009      FORMAT ( f  T,TCOEFF.  OF  Z**( ',12, ')=' ) 
READ( 5,*)  C6(  I ) 

2008  CONTINUE 

DO  112  L=l,200 

VP( L)=SIN( L*PI/3)+SIN( L*PI/4)+SIN( L*PI/5) 
+      +  SIN( L*PI/7)+SIN( L*PI/2)+SIN(  L*PI/6) 
112    CONTINUE 

DO  243  L=l,200 
VC(  L)=100.  0 
243    CONTINUE 

DO  124  L=1,N 
Y(L)=0.  6 
U(LJ=VP(L) 
124    CONTINUE 
Y(N1}=0.  0 
DO  2010  1=1, N 
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™,«      „      Y(N1)=Y(N1)-C6(I)*Y(I)+C5(I)*U(I) 

2010  CONTINUE 

DO   2011    1  =  1, N 
J9=N1-I 
L5=N6-I 

U3( J9,M)=C6( I) 
U3    L5,M)=C5( I 

2011  CONTINUE 

DO  113  1=1, N 
L1=N+I 
L2=N1-I 


Fl( I,M)=0.  0 
Fl(Ll,M)=VP(L2 


113    CONTINUE'  ' 

DO  121  1=1. N2 

DO  122  K=1,N2 

IF(  I.  EQ.  K)  THEN 
P2(  I,K)=30. 0 
ELSE 

P2( I,K)=0.  0 
END  IF 

122  CONTINUE 
121    CONTINUE 

C      CALL  FAMILY(N2,N2, ' P2 ' ) 
C      CALL  RESULT(N2,N2, I.K.P2) 

CALL  COPYM(N2/N2/ I,L,AB,P2) 

DO  123  1=1, N2 
Ul( I,M)=0.  0 

123  CONTINUE 

CALL  TRANS(N2.M, I,K,F1,FT) 
CALL  FRTCMS(  TCLRSCRN  T) 
WRITE(6.702) 

702  FORMATC   T,fENTER  Q(Z)  IN  DESCENDING  ORDER  OF  Z') 
WRITE(6.703) 

703  FORMAT(4  r, 'HIGHEST  COEFF.  MUST  BE  1.  ' ) 
WRITE(_6.704) 

704  FORMAT ( f  T , f *REMARK:  Q( Z)  SHOULD  BE  STABLE  POLYNOMIAL') 

DO  157  i=l,N 

WRITE(6.705)N-I 

705  FORMAT (  '  T/COEFF.  OF  Z**(  '  ,  12,  '  )  =  '  ) 
READ(5,*)  C3(  I) 

157       CONTINUE 

CALL  FRTCMS( ' CLRSCRN  ') 
WRITEf6.706) 

706  FORMAT ( '  T, 'ENTER  PSTAR( Z)  IN  ASCENDING  ORDER 

OF  Z'  ) 
WRITE( 6.707) 

707  FORMAT(*  T, 'HIGHEST  COEFF.  MUST  BE  1. ' ) 
WRITEr6.708) 

708  FORMAT (  *  T ,  ' * REMARK- 1:  PSTAR(  Z )  IS  THE  DENOMINATOR 
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709 
710 

711 
158 

241 


116 

115 


POLYNOMIAL1 ) 
WRITE(6.709) 
FORMATS '  T,'OF  MODEL1) 
WRITE(  5f 7}0) 


FORMAT (  '  '  ' *REMARK-2:  PSTAR(Z)  SHOULD  BE  STABLE 

POLYNOMIAL1  ) 
DO  158  1=1, N 

OF  Z**( ' , 12,  * )  =  ' ) 


WRITE(6.711)N-I 

FORMAT( '  T . 'COEFF. 

READ(5,*)  C4(  I) 
CONTINUE 
DO  241  1=1,20 

V(  I)=VP(l) 
CONTINUE 
DO  114  J=N1,L8 

CALL  CALCUL( N2 , M, N2 , I , L . K, VI, P2 , Fl ) 

CALL  CALCUL(M,M/N2/ I , L, K, V2 , FT, VI ) 

TSP=1. +V2CM.M) 

CALL  CALCUL(M/M,N2,I,L,K/V4,FT,U1) 

TSF=Y( J)-V4(M,M) 

TSN=TSF/TSP 

CALL  SCALAR(N2,M. I , K, TSN, VI, V5 ) 

CALL  SUM(N2,M. I , K, Ul, V5,U2 ) 

CALL  CALCUL(N2yM,N2, I , L, K, P3 , P2 , Fl ) 

CALL  CALCUL(M, M, N2 , I , L, K, P4, FT, P3 ) 

TSM=l./(  1  +  P4(M,M)  ) 

CALL  CALCUL(M.N2.N2, I , L, K, P5, FT, P2 ) 

CALL  CALCUL( N2 , N2 , M . I , L , K . P6 , Fl . P5 ) 

CALL  CALCUL  N2 , N2 , N2 , I , L, K, P7 , P2 , P6 ) 

CALL  SCALAR  N2,N2, I, L,TSM,P7/P8 ) 

CALL  SUB(N2,N2/I/L,P2,P8,P9) 

INITIALIZE 

DO  115  1=1. N2 

DO  116  K=1,N2 
P2( I,K)=0. 0 

CONTINUE 


CONTINUE 
CALL  COPYM(N2,N2 
IF( J. EQ. 10)  THEN 
CALL  COPYM(N2/LN2 
ELSE  IF?J.  EQ.  20) 
CALL  COPYM(N2.N2 
ELSE  IF? J. EQ. 30) 
CALL  COPYM(N2,N2 
ELSE  IF(J.  EQ.  40) 
CALL  COPYM(N2,N2 
ELSE  IF(  J.  EQ.  50) 
CALL  COPYM(N2,N2 
ELSE  IF( J. EQ.  60) 
CALL  COPYM(N2,N2 


I,L,P2,P9) 

I,L,P2,AB) 

THEN 

I,L,P2,AB) 

THEN 

I,L,P2,AB) 

THEN 

I,L,P2,AB) 

THEN 

I,L,P2,AB) 

THEN 

I,L,P2,AB) 
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ELSE  IF(J.EQ.  70]  THEN 

CALL  COPYM(N2/N2/I/L/P2,AB) 

END  IF 

IF( J. LE. M7)  GO  TO  9591 

C5(l)=3. 
DO  2031  1  =  1, N 

J9=N1-I 

L5=N6-I 

U3( J9,M)=C6(  I) 

U3( L5,M  =C5  I) 
2031   CONTINUE 
9591      DO  117  1=1, N2 

Ul( I,M)=0.  0 
117       CONTINUE 

CALL  C0PYM(N2.M,I,L,U1,U2) 
C         CALL  FAMILY(N2,N2,rP2Tj 
C         CALL  RESULT  N2,N2, I, K,P2) 

DO  148  1  =  1, N 
U( I)=V( i) 
148       CONTINUE 

DO  211  1=1, N 

D( I)=U1( I,M) 

211  CONTINUE 

DO  212  1=1, N 
L4=N+ I 
B(  I)=U1(L4,M) 

212  CONTINUE 

DO  7011  I=1,N 

M9=N-I+1 

Ql( I)=C3(M9) 
7011      CONTINUE 

Nl)=l.  0 

N1]=0. 0 

4001  1=1, N 

Q3( I)=D( t)-C4( I) 
4001      CONTINUE 

DO   7007    I=1,N 

M8=N-I+1 

nrsnn  _XT§?  (  I)  =Q3(  M8  ) 

7007      CONTINUE 
DEG1=N+1 
DEG2=2*N 
DEG3=DEG2+2 
DO  6001  I=1,DEG3 
R(I)=0. 0 
DO  6002  K=1,DEG2 

IF(  (  (.I-K).  LT.  0.  0).  OR.  (  (  I-K).  GT.  DEG1 )  ) 

GO  TO  6002 
IF((  I.GT.  l).OR.  (K.  GT.  1)}  GO  TO  6503 
6503  R(  I-1)=R( I-1)+Q1(K-1)*Q2(  I-(K-l) ) 


n 
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6002         CONTINUE 
6001      CONTINUE 

R(1)=Q1(  1)*Q2( 1) 
DO  9003  I=1*N6 
K9=N6-I+1 
F(  I)=R(K9) 
9003      CONTINUE 
N4=2*N+1 
M1=N+1 

DO  41  J4=1/M1 
DO  51  1=1, N4 
K=I-J4 


EQ. 0.)    THEN 


A(I,34) 


HW-0 

)  IF 

IF  ((K.  GT.  0.  ).  AND.  (K.  LT.M1))  THEN 
Xtl/J4)=D(K) 

END  IF 
51        CONTINUE 
41     CONTINUE 

DO  61  J4=1,N 

DO  71  1=1. N4 

A( I/N+l+J4)=0. 

K=I-J4-1 

IF  ((K.  GT.  0.  ).  AND.  (K.  LT.  Ml)  )  THEN 

A(.I/N+1+J4)=B(K) 
END  IF 
71        CONTINUE 
61     CONTINUE 

DO  10  C0LUM=1/N4 

SUM2=0. 

DO  70  J4=C0LUM,N4 

SUM2=SUM2  +  (A(  J4/C0LUM) )**2 
70     CONTINUE 

SIGMA=SQRT(SUM2) 

ASIGMA=ABS( SIGMA) 
IF  ( ASIGMA.  LT.  0.  1)  GO  TO  6666 
C  WRITE(6.714) 

C14  FORMAT ( '  r7'MATRIX  A  IS  SINGULAR") 

C         END  IF 

DO  80  J4=COLUM,N4 
G(  J4J=0. 
80     CONTINUE 

G(COLUM)=SIGMA 

SUM2=0. 

DO  90  J4=C0LUM,N4 

SUM2=SUM2  +  (A(  J4,C0LUM)-G(  J4) )**2 
90     CONTINUE 

DO  100  J4=C0LUM7N4 
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DO  110  K=C0LUM,N4 

IF  (ABS( SUM2 ) . LT. 0. 00001 )  THEN 

TEMP=0. 
ELSE 

TEMP=-2*{ A( J4, COLUM) -G( J4) ) *(  A( K, COLUM) 
+  -G(K))/SUM2 

END  IF 
IF  (  J4. EQ. K)  THEN 

TEMP=I.  +  TEMP 
END  IF 

E( J4,K)=TEMP 
110    CONTINUE 
100    CONTINUE 

DO  131  J4=COLUM,N4 
DO  130  K=C0LUM,N4 
TEMP=0. 
DO  140  L=C0LUM,N4 

TEMP=TEMP+H(  J4 , L ) * A(  L , K ) 
140         CONTINUE 

C[J4,K)=TEMP 

130  CONTINUE 

131  CONTINUE 

DO  150  J4=C0LUM/N4 
TEMP=0. 
DO  160  L=C0LUM.N4 

TEMP=TEMP+H( J4, L ) *F(  L) 
160      CONTINUE 

T(_  J4 )  =TEMP 
150    CONTINUE 

DO  170  J4=COLUM,N4 
DO  180  K=C0LUM,N4 
A( J4,K)=C(  J4,K) 
180      CONTINUE 

F(_J4J=T(  J4) 
170    CONTINUE 
10     CONTINUE 

X( N4 ) =F( N4 ) /A( N4 , N4 ) 
I=N4-1 
812      SUM2=0. 
M2=I+1 
DO  30  J4=M2.N4 

SUM2=SUM2  +  A(  I, J4)*X( J4) 
30       CONTINUE 


X(I)=(F(  I)-SUM2)/A(  I.I] 

IF  (ABS(X(  I)  ).  LT.  0.00001)  THEN 

X( I)=0. 
END  I! 

1  =  1-1 

IF(  I.  i 

6666     HI  J)^ 


END  IF 
1  =  1-1 

IF( I. GE. 1)  GO  TO  812 
>=X(1) 
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S2    I)=X(I 

141  CONTINUE 

DO    142    I=N3,N4 
N5=I-N1 
S1(N5)=X(  I) 

142  CONTINUE 
UP(  J)=0.  0 

DO  5001  I  =  1,N 
L7=I+1 

UP( J)=UP[J)+(S2(L7)-C3(  I))*U(  J-I)+S1(  I)*Y(  J-I) 
+  +C3( I)^V( J-I ) 

5001     CONTINUE 

U(J)=( UP( J)+V( J) )/( 1-S2(  1) ) 
J5=J+1 
Y(J5}=0.  0 
DO  3003  I=1#N 
M6=N1-I 
Y( J5 )=Y( J5 ) -C6( I ) *Y( J5-M6 ) +C5(  I ) *U(  J5-M6 ) 

3003  CONTINUE 

DO  189  1=1, Nl 

YM( 11=0. 0 
189     CONTINUE 

YM(  J5}=0.  0 
DO  2508  I=1,N 

M8=N-I+1 

Bl(  I)=C5(M8) 
2508    CONTINUE 

DO  3004  1  =  1,  N 

YM( J5)=YM( J5)-C4( I )*YM( J5-I )+Bl( I)*V( J5-I) 

3004  CONTINUE 

MK(  J5)=-Y( J5) 
DO  3005  1=1. N 

MK( J5)=MK( J5)+B( I )*V( J5-I )-C4(  I)*Y(  J5-I) 

3005  CONTINUE 

MU(  J5)=ABS(MK(_J5)) 
IF(MU(J5).  LT.  20.  0)  THEN 

V(J5)=^C(J5) 
ELSE 

V(J5)=VC(J5)+VP(J5) 

END  IF 

DO  120  I=1,N2 
Fl( I,M)=0. 0 
120       CONTINUE 

DO  118  I=1,N 
L1=N+I 
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I2AM,I,K.F1,FT) 

n2,m,  *uiT  ) 

N2,M, I,K,U1) 


J2=J+1-I 

Fl(  I.M_)=-Y(  J2) 

F1(L1,M)=U(  J2 

118  CONTINUE 
DO  119  1=1, N2 

FT(M/I)=0.  0 

119  CONTINUE 
CALL  TRANS(N2 
CALL  FAMILY( 
CALL  RESULT 
DO  3001  1=1, N2 

U4(  I,M)=0.  0 

3001  CONTINUE 

CALL  SUBCN2,N2, I,L,U3,U1,U4) 

C7( J)=0. 0 

DO  3002  1=1 ,N2 

C7( J)=C7( J)+U4( I,M)*U4(  I,M) 

3002  CONTINUE 

C8( J)=SQRT(C7( J) ) 
114       CONTINUE 

DO  8888  I=17N1 
C8( I )=0. 
8888      CONTINUE 
C         DO  8111  1=1,100 
C         WRITE  j 8,8169)  C8(  I) 
Clll      CONTINUE 
C169      FORMATCIOX.FI?.  6) 

WRITE(8,lllir  SYS.  OUTPUT' ,  'MOD.  OUTPUT 
1111      FORMAT( 13X,A11,6X,A11) 

DO  976  1=1,100 
C         WRITE  (6,8002)  Y(I)/YM(I 

WRITE  (8,8002  Y(.  I )  ,  YM(  I 
8002  FORMAT] 9X,F15. 6,6X,F15. 6 
976       CONTINUE 

DO  4441  1=1, L8 
C  WRITE( 8,4445)1 

C  WRITE  8,4442)  Hl( I 

C  WRITE  8,4442)  H2( I 

C  WRITE  8,4442)  H3( I 

C  WRITE  8,4442)  H4(  I 

C  WRITE( 8,4442)  H5( I 

4445      FORMAT(T  T, 12) 
4442      FORMAT? 10X,F15.  6) 
4441      CONTINUE 

DO  1025  I=1,L8 

SE  i;2)=YM(  I) 

KD  11=1 

1025      CONTINUE 

r   ****************************************** 
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C   *  THIS  PART  PLOTS  THE  NECESSARY   DATA     * 
C   *  THAT  PROGRAM  CALCULATES  * 

Q    ****************************************** 

CALL  VPLOT(KD,SE,L8,2, 'TIME 
CALL  Vl?LOT(  KD, C8,L8,1, [TIME 


CALL  V2PLOT  KD,MU, L8.1. 'TIME 
CALL  V3PLOT  KD,V,L8,1, 'TIME 
CALL  V4PL0T  KD,U,L8,1, 'TIME   ' 
CALL  V4PL0T  KD,H1,L8, 1, ; 


i 
'  i 

i 
i 

t 
'  i        i 

'  i        i 
'  »        i 

'  !  » 


C 

C  CALL 
C  CALL 
C         CALL 

c      call  v4plot(kd;h5;l8;i; 

STOP 
END 

Q    ****************************************** 

C   *  THIS  SUBROUTINE  PRINTS  THE  ARRAYS  AS  AN* 
C   *  MATRIX  FORM   -  * 

Q  ****************************************** 

SUBROUTINE  RESULT( H, O, P, S, T) 

INTEGER  H  ,  O  ,  P  ,  S 

REAL  T  (  H  ,  0  ) 

DO  80   P  =  1  ,  H 

WRITE( 8, 70)(T(  P, S) , S=l, O) 
C         PRINT  70,  (T(_P,  S_),  3=1.0) 
70        FORMAT(  T6\T2,9X,10(3X,F12.  6)) 
80      CONTINUE 

RETURN 

END 
r*      ****************************************** 

C   *  THIS  SUB.  CALCULATES  THE  MULTIPLICATION* 
C   *  TWO  MATRICES.  * 

Q    ****************************************** 

SUBROUTINE  CALCUL( U, V, Y, Z, X, ZX, ZT, W, Q) 
INTEGER  U,V,Y, ZVX, ZX 
REAL  ZT(U,V  ,W(U,Y),Q(Y,V) 
DO  93  Z=1,U 
DO  92  X=1,V 

ZT(Z,X)  =  0. 0 
DO  91  ZX=1,Y 

ZT(Z  X)=ZT(Z,X)+W(Z,ZX)*Q(ZX,X) 

91  CONTINUE 

92  CONTINUE 

93  CONTINUE 
RETURN 
END 

Q    ****************************************** 

C   *  THIS  SUBROUTINE  WRITES  THE  MATRIX  NAME  * 
C   *  AND  ROW. COLUMN  NUMBERS.  * 

Q    ****************************************** 


J 
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subrout i ne  fam i ly( row , column , name ) 
integer  row, column 
character* 2  name 
write(8,94)  ^ matrix'  name 
Matrix1  .name 


C       PRINT  94 


94  FORMATfTi, '0f,T2, i0X; T12 . A6, T19, IX, T20, A2 ) 
WRITE(8,95)  *  ROW  NUMBER.  , ROW 

C       PRINT  95/ROW  NUMBER.  'ROW 

95  FORMATf Tl,  0' ,T2, 10X, T12 , All, T23 , 8X, T31, 12) 
WRITE(8,96)  'COLUMN  NUMBER. T, COLUMN 

C       PRINT  96.rCOLUMN  NUMBER. r. COLUMN 

96  FORMAT( Tl , r0T, T2 , 10X, T12 , A14, T26, 5X, T31, 12 ) 

RETURN 

END 
p   ******************************************** 

C   *  THIS  SUBROUTINE  CALCULATES  THE  SUMMATION  * 
C   *  OF  THE  SAME  COLUMN  ELEMENTS.  * 

q      ******************************************** 

SUBROUTINE  ALI( II , 12 , HI, H2 , Gl, W2 ) 
INTEGER  II, 12, HI, H2 
REAL  W2(l, 12), Gl( 11.12) 
DO  41  H2=l, 12 

DO  42  Hl=l, II 
W2( 1,H2)=W2( 1/H2)+ABS(G1(H1,H2)) 

42  CONTINUE 
41      CONTINUE 

RETURN 
END 

Q         ****************************************** 

C   *  THIS  SUBROUTINE  MULTIPLIES  THE  MATRIX   * 
C   *  BY  SCALAR  NUMBER.  * 

Q  ****************************************** 

SUBROUTINE  SCALAR( Kl , K2 , Gil , G2 , PI , G3 , G4) 
INTEGER  K1,K2,G11,G2 
REAL  G3(K1,K2) . G4(  Kl, K2 ) , PI 
DO  43  G2=1,K2 

DO  44  G11=1,K1 
G4(G11,G2)=G3(G11,G2)*P1 
44      CONTINUE 

43  CONTINUE 

RETURN 

END 
C   ****************************************** 

C   *  THIS  SUBROUTINE  CALCULATES  THE  SUM  OF   * 
C   *  THE  TWO  MATRICES.  * 

Q    ****************************************** 

SUBROUTINE  SUM(K3yK4, 15, K5 ,X1,X2 ,X3 ) 
INTEGER  K3,K4, I5.K5 

REAL  Xl(  K3  ,  K4_)  ,  X2(  K3  ,  K4)  ,  X3(  K3  ,  K4) 
DO  45  K5=1,K4 


101 


DO  46  15=1, K3 

X3( I5,K5)=X1( I5/K5)+X2(  I5,K5) 
46      CONTINUE 

45  CONTINUE 
RETURN 
END 

Q  ****************************************** 

C   *  THIS  SUBROUTINE  COPPIES  THE  MATRICES.   * 

Q  ****************************************** 

SUBROUTINE  COPYM(  K3 , K4, 15 , K5 , XI , X2 ) 

INTEGER  K3,K4, I5.K5 

REAL  X1(K3,K4J,X2(K3,K4) 
DO  45  K5=1,K4 

DO  46  15=1. K3 

Xl( I5,K5)=X2(  15, K5) 

46  CONTINUE 

45  CONTINUE 
RETURN 
END 

p   ****************************************** 

C   *  THIS  SUBROUTINE  FINDS   THE  TRANSPOSE-   * 
C   *  OF  THE  MATRIX.  * 

Q    ****************************************** 

SUBROUTINE  TRANS(K3/K4/ I5/K5/X1/X2) 

INTEGER  K3 , K4 , I 5 . K5 

REAL  XI (  K3  ,  K4_) ,  X2  (  K4 ,  K3  ) 
DO  45  K5=1,K4 

DO  46  15=1. K3 
X2(K5/I5)=X1(  15^5) 

46  CONTINUE 

45  CONTINUE 
RETURN 
END 

Q    ********************************************** 

C   *  THIS  SUBROUTINE  CALCULATES  THE  SUBTRACTION.* 

C   *  OF  THE  MATRICES.  * 

C   ********************************************** 

SUBROUT INE  SUB(  K3 . K4 , 1 5 , K5 , XI , X2 , X3 ) 
INTEGER  K3  ^4.15.^5 

REAL  XI (  K3  ,  K4_)  ,  X2  (  K3  ,  K4 )  ,  X3  (  K3  ,  K4 ) 
DO  45  K5=1,K4 

DO  46  15=1, K3 

X3(  I5,K5)=X1( I5,K5)-X2( 15, K5) 

46  CONTINUE 
45      CONTINUE 

RETURN 
END 
C   ****************************************** 

C   *  THIS   SUB.  PLOTS   THE   ARRAYS   USING    * 
C   *  THE  DISSPLA  PROGRAM.  * 
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c\        ****************************************** 

C 

SUBROUTINE  VPLOT( T,U,M, N, LABELX, LABELY) 

C************* 

REAL  U(  800,4) ,T( 800) ,Y( 800) .YX(4),YN( 4) 
INTEGER  PMAX, PMIN,M,N, LABELX,  LABELY 
C 

CALL  AMAX(T,M,TMAX,PMAX) 
CALL  AMIN{T,M,TMIN,PMIN) 
DO  20  J  =  1,N 

DO  10  I  =  1,M 


Y(I)  =u(l,J) 
ITINUE 


io     cont: 

CALL  AMAX(  Y , M , YMAX , PMAX ) 
YX(  J)  =  YMAX 
CALL  AMIN(Y,M,YMIN,PMIN) 
YN( J)  =YMIN 
WRITE(6,11)  YMAX,YMIN 
20    CONTINUE 

WRITE( 6/13)(YX( I), 1=1,4) 
WRITE( 6. 13J ( YN(  IJ  . 1  =  1, 4J 
CALL  AmAX(YX,N, YMAX, PMAX) 
CALL  AMIN[YN,N.YMIN,PMIN) 
WRITE(6, 11}  YMAX,YMIN 
CALL  TEK  618 
CALL  BL0WUP(1. 2) 
CALL  PAGE(  11.  ,8.  5) 
:      CALL  NOBRDR 

CALL  AREA2D(9.  ,6.  ) 

CALL  XNAME( LABELX, 6 ) 

CALL  YNAME( LABELY, 6 ) 

CALL  HEADIN  ( TPLANT  AND  MODEL  OUTPUTS  $',100,1.,1) 

CALL  CROSS 

CALL  GRAF( TMIN, ' SCALE ' , TMAX, YMIN, ' SCALE ' , YMAX) 

CALL  SPLINE 

DO  40  J  =1,N 

DO  30  I  =1,M 

Y(I)  =U(l,J.) 
30       CONTINUE 

IF( J. EQ. 2 )  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 
40    CONTINUE 

CALL  ENDPL(O) 
:      CALL  DONEPL 

11   FORMATS1  ;  , 10X,2G12.  5/) 
13   FORMAT (  "  '  ,10X,4G12.  5/) 
RETURN 
END 

SUBROUTINE  V1PL0T( T, U, M, N, LABELX, LABELY) 
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REAL  U( 800,4) ,T( 800) , Y( 800 ) , YX( 4), YN( 4) 
INTEGER  PMAX,PMIN,M,N, LABELX, LABELY 
C 

CALL  AMAX(T,M,TMAX,PMAX) 
CALL  AMIN(T,M,TMIN,PMIN) 
DO  20  J  =  1,N 

DO  10  I  =1,M 

Y(I)  =u(l,J) 
10       CONTINUE 

CALL  AMAX(Y,M,YMAX,PMAX) 
YX( J )    =   YMAa 
CALL  AMIN(Y,M,YMIN,PMIN) 
YN( J)  =YMIN 
WRITE(6,11)  YMAX,YMIN 
20    CONTINUE 

.  WRITE(6/13)(YX(  I)/I  =  l/4) 
WRITE(6.13J(YN( I) .1=1,4 j 
CALL  AMAX(YX/N/YMAX/PMAX) 
CALL  AMIN(YN,N.YMIN,PMIN) 
WRITE(6. 11)  YMAX,YMIN 
CALL  TEK  618 
CALL  BLOWUP(l.  2) 
CALL  PAGE( 11.  ,8.  5) 
C      CALL  NOBRDR 

CALL  AREA2D(9.  ,6.  ) 

CALL  XNAME( LABELX, 6 ) 

CALL  YNAME(  LABELY. 6J 

CALL  HEADIN  ( TPARAMETER  ERROR  $' ,100,1. ,1) 

CALL  CROSS 

CALL  GRAF(TMIN, 'SCALE' , TMAX , YM I N , 'SCALE' ,YMAX) 

CALL  LINEAR 

DO  40  J  =1,N 

DO  30  I  =  1,M 

Y(IJ  =U(I,J) 
30       CONTINUE 

IF{ J. EQ.  2  )  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 
40    CONTINUE 

CALL  ENDPL(O) 
C      CALL  DONEPL 

11   FORMAT(^  1 ,10X,2G12.  5/) 
13   FORMATJ'  ' ,10X,4G12.  5/) 
RETURN 
END 
C 

SUBROUTINE  V2PLOT(  T, U,M,N, LABELX, LABELY) 

Q************* 

REAL  U(  800.4)  ,T(  800) ,Y(  800). YX( 4), YN( 4) 
INTEGER  PMAX, PMIN, M, N, LABELX, LABELY 
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CALL  AMAX(T,M,TIYIAX,PMAX) 
CALL  AMIN(T,M/TMIN/PMIN) 
DO  20  J  =  1,N 

DO  10  I  =1,M 

Y(I)  =U(I,J) 
10       CONTINUE 

CALL  AMAX( Y , M , YMAX , PMAX ) 

YX(  <J )    =  YMAa 

CALL  AMINtY/M/YMIN/PMIN) 

YN( J)  =YMIN 

WRITE(6,11)  YMAX,YMIN 
20    CONTINUE 

WRITE( 6#13)(YX(I)# 1=1,4) 
WRITE(6,  13  W  YN  IJ  ,  1=1, 4  J 
CALL  AmAx(YX,N, YMAX, PMAX) 
CALL  AMIN(YN.N,YMIN,PMIN 
WRITE(6.11]  YMAX,YMIN 
CALL  TEK  618 
CALL  BLOWUP(l.  2) 
CALL  PAGE( 11.  ,8.  5) 
:      CALL  NOBRDR 

CALL  AREA2D(9.  Jo.  ) 

CALL  XNAME( LABELX, 6) 

CALL  YNAME( LABELY, 61 

CALL  HEADIN  ( TOUTPUT  PREDICTION  ERROR  $' ,100,1. ,1) 

CALL  CROSS 

CALL  GRAF(TMIN, 'SCALE' , TMAX , YM I N , 'SCALE' ,YMAX) 

CALL  SPLINE 

DO  40  J  =1,N 

DO  30  I  =  1,M 
Y(IJ  =U(I,J) 
30       CONTINUE 

IF(J.  E0.2)  CALL  CHNDOT 

CALL  CURVE  (T,Y,M,0) 
40    CONTINUE 

CALL  ENDPL(O) 


CALL  DONEPL 

FORMAT* ^  1 ,10X,2G12.  5/) 
13   FORMAT)'  '  ,10X,4G12.  5/J 
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RETURN 
END 


SUBROUTINE  V3PLOT( T,U,M,N, LABELX, LABELY) 

REAL  U( 800.4 ),T( 800) , Y(  800) , YX( 4) , YN( 4) 
INTEGER  PMAX, PMIN,M,N, LABELX, LABELY 
C 

CALL  AMAX(T,M,TMAX,PMAX) 
CALL  AMIN(T/M/TMIN/PMIN) 
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DO  20  J  =1,N 

DO  10  I  =1,M 
Y(I)  =u(l,J) 
10       CONTINUE 

CALL  AMAX(Y.M,YMAX,PMAX) 

YX(  J )    =  YMAX 

CALL  AMIN(Y,M,YMIN,PMIN) 

YN( J)  =YMIN 

WRITE(6,11)  YMAX#YMIN 
20    CONTINUE 


CALL  AMIN(YN,N.YMIN,PMIN) 
WRITE(6. 11)  YMAX,YMIN 
CALL  TEK  618 
CALL  BLOWUPfl.  2) 
CALL  PAGEfll:  ,8.  5) 
C      CALL  NOBRDR 

CALL  AREA2D(9.  ,6.  ) 

CALL  XNAME( LABELX, 6) 

CALL  YNAME  LABELY, 6J 

CALL  HEADIN  ( T EXTERNAL  INPUT  $',100,1. ,1) 

CALL  CROSS 

CALL  GRAF(TMIN, 'SCALE' , TMAX, YMIN, 'SCALE' ,YMAX) 

CALL  SPLINE 

DO  40  J  =1,N 

DO  30  I  =1,M 
yTi)   =U(I,J) 
30       CONTINUE 

IF(J.  EQ.2)  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 
40    CONTINUE 

CALL  ENDPL(O) 
C      CALL  DONEPL 

11   FORMAT^  7 ,10X,2G12.  5/) 
13   FORMAT}'  ' ,10X,4G12.  5/j 
RETURN 
END 
SUBROUTINE  V4PL0T( T, U, M, N, LABELX, LABELY) 

REAL  U( 800 . 4 ) , T( 800 ) , Y[ 800) . YX( 4 ) , YN( 4 ) 
INTEGER  PMAX, PMIN, M, N, LABELX, LABELY 
C 

CALL  AMAX(T,M, TMAX, PMAX) 
CALL  AMIN[T,M,TMIN,PMIN) 
DO  20  J  =1,N 

DO  10  I  =1.M 
YflJ  =U(I,J) 
10       CONTINUE 
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CALL  AMAX(Y,M, YMAX, PMAX) 
YX( J )  =  YMAX 
CALL  AMINCY/M/YMIN/PiyiIN) 
YN(  J)  =YMIN 
WRITE(6,11)  YMAX,YMIN 
20    CONTINUE 


11 
13 


CALL  AMIN(YN,N,YMIN,PMIN 
WRITE(6.1l]  YMAX,YMIN 
CALL  TEK  618 
CALL  BLOWUP(l. 2) 
CALL  PAGEfll.  ,8.  5) 
C  ■     CALL  NOBRDR 

CALL  AREA2D(9.  ,6.  ) 

CALL  XNAME( LABELX, 6 ) 

CALL  YNAME( LABELY, 6 ) 

CALL  HEADIN  ( TINPUT  TO  THE  PLANT  $',100,1. ,1) 

CALL  CROSS 

CALL  GRAF(TMIN, 'SCALE' , TMAX, YMIN, 'SCALE' ,YMAX) 

CALL  SPLINE 

DO  40  J  =1,N 

DO  30  I  =1,M 

Y(I)  =U(I,J) 
30       CONTINUE 

IF(J. EQ. 2)  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 
40    CONTINUE 

CALL  ENDPL(O) 
C      CALL  DONEPL 

FORMAT(^  7 ,10X,2G12.  5/) 
FORMAT (  '  '  ,10X,4G12.  5/) 
RETURN 
END 


SUBROUTINE  AMAX( Y,N, YMAX, PMAX) 

Q  *****************************************  it.-k  -k-k  -k  -k 

C         *   SUBROUTINE  TO  COMPUTE  THE  LARGE  ELELEMENT 
C         *   IN  A  GIVEN  ROW  OR  COLUMN  IN  ARRAY  TAT. 

c 

C         ***  VARIABLE  DECLARATION  *** 
C         REAL  AMAX 

INTEGER  PMAX 

DIMENSION  Y(N) 


YMAX=Y(  1) 


DO  5100  1=2, N 
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IF(Y(  I).  LE.  YMAX)  GO  TO  5000 
YMAX=Y( I ) 
PMAX=I 
5000         CONTINUE 
5100      CONTINUE 
C         MAXR0W=P0S 
RETURN 
END 
C 
C 

SUBROUTINE  AMIN( Y,N, YMIN, PMIN) 
C 

C         *   SUBROUTINE  TO  COMPUTE  THE  SMALL  ELEMENT 
C         *   IN  A  GIVEN  ROW  OR  COLUMN  IN  ARRAY  TAT. 

Q  ******************************************** 

c 

C         ***  VARIABLE  DECLARATION  *** 
C         REAL  AM IN 

INTEGER  PMIN 
DIMENSION  Y(N) 
C 

YMIN=YC1) 

DO  6100  1=2, N 

IF(Yfl).GE. YMIN)  GO  TO  6000 
YMIN=Y( I ) 
PMIN=I 
6000  CONTINUE 

6100      CONTINUE 
RETURN 
END 
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APPENDIX  D 
COMPUTER  PROGRAM  RCIOP 

P    ****************************************** 

C  *  THESIS  PROGRAM  * 

C  *  INDIRECT  ADAPTIVE  CONTROL  * 

C  *  PARAMETER  ESTIMATION  * 

C  *  USING  PROJECTION  ALGORITHM  * 

Q  ****************************************** 

C        DEFINITION  OF  VARIABLES: 

INTEGER  N1,N,L,K,L1. I,M,L2, J/N3/N2/ Jl, J2 , J3/M1/M2 
INTEGER  N4, N5 , N6 . COLUM, J5 , M5 , J7 , J8 , J9 , L5 , M6 , L6 , L7 
INTEGER  DEG1 , DEG2 , DEG3 . K9 , M7 , L8 , M9 , M8 , M3 , L4 
REAL  V1(4,1)\V2£  1,1),V3(  4.  1).V4(  l.i),V5(4,l) 
REAL  U2(4,1),Y(  800)  .U(  800)  .Fl(  4,  1}  .FT(  1,4)  .P2(  4.4) 
REAL  AB(4,4J  ,P3(  4,4)  ,P4(4,  1 J i,  81(20)  ,U1(4.  1 ) ,  R  20) 
REAL  TSP,TSF,TSN.PI .V( 800), YM( 800) ,MU( 800), VP  800) 
REAL  A( 20,20) . H( 20, 20) , B( 20) , C( 20, 20} . X( 20) 
REAL  D£10KT(20)  .F(20),G}20J  .Sl(  10),S2(  10),C3( 10) 
REAL  U3(  4 '  1 )  ,  C5(  10}  ,  C6(  10}  .  MK(  800)  ,  U4(  4,1,  C7  800  ) 
"D(806J,C8(806),C4(  l65,Q3(50) 
, SIGMA, ASIGMA,DY.UP( 800}, Q2< 
,se(806,4),v6(  l,l),ql(50) 


REAL  VC  86oj ,KD 

REAL  SUM2,TEMP, SIGMA, ASlGMA! DY.UP( 800);Q2( 50) 
REAL  CON1.CON2 
PI=3. 141592654 
DO  20  I  =1,200 
DO  8   J  =1,4 


SE(  I, J)=0. 0 
8       CONTINUE 
20    CONTINUE 
M=l 

WRITE(_6,701) 
701    FORMAT(  ' lr  ,  f ENTER  THE  ORDER  OF  THE  SYSTEM') 
READ(5,*)  N 
N1=N+1 
N2=2*N 
N3=N+2 
N6=2*N+1 
L8=10*N2 
WRITE( 6.2301) 

2301  FORMAT ( '  T,TENTER  CONSTANT  C  (PROJECTION  ALGORITM) ' ) 
READ (5,*)  CON1 

WRITER.  2302) 

2302  FORMAT(  f  T,TENTER  CONSTANT  A  (PROJECTION  ALGORITHM)') 
READ(5,*)  CON2 

WRITE( 6.2002) 
2002   FORMAT( flT,T ENTER  PLANT  NUMERATOR  POLYNOMIAL 
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COEFFICIENTS' ) 
WRITE(_6.2003) 

2003  FORMAT (  '  T.TIN  ASCENDING  ORDER  OF  Z') 
DO  2001  I  =  1,N 

J7=I-1 

WRITE( 6.2004)  J7 

2004  FORMAT ( f  T,rCOEFF.  OF  Z**l  ' , 12, ' )  =  ' ) 
READ(5,*)  C5(  I) 

2001   CONTINUE 

WRITE(  6.2005) 

2005  FORMAT( flT,T ENTER  PLANT  DENOMINATOR  POLYNOMIAL 

COEFFICIENTS7) 
WRITE(6.2006) 

2006  FORMATC   T,rIN  ASCENDING  ORDER  OF  Z ' ) 
WRITE(6.2007) 

2007  FORMATC    ,  HIGHEST  COEFF.  SHOULD  BE  1.  ') 
DO  2008  1=1, N 

J8=I-1 

WRITE( 6.200?)  J8 

2009  FORMAT (  i    T,rCOEFF.  OF  Z**(  *  ,  12 ,  '  )  =  '  ) 
READ(5>)  C6(  I) 

2008  CONTINUE 

DO  112  L=l,200 

VP( L)=SIN( L*PI/3)+SIN( L*PI/4)+SIN(  L*PI/5) 
nio   +^^T  +  sin(l*pi/6)  +  sin(l*pi/7)  +  SIN  L*PI/2 

112  CONTINUE 

DO  243  L=1.200 

VC( L)=100. 0 
243    CONTINUE 

DO  124  L=1.N 

Y(L)=0.  6 

U(LJ=VP(L) 
124    CONTINUE 
Y(N1)=0.  0 
DO  2010  I  =  1,N 

Y(N1)=Y{N1)-C6(  I)*Y(  I)+C5(  I )  *U(  I ) 

2010  CONTINUE 

DO  2011  I=1,N 
J9=N1-I 
L5=N6-I 

U3( J9,M)=C6(  I) 
U3(L5,M)=C5(  I) 

2011  CONTINUE 

DO  113  I  =  1,N 
L1=N+I 
L2=N1-I 
Fl(  I.M)=0.  0 
Fl  Ll,M)=VP(L2) 

113  CONTINUE 

DO  121  I=1,N2 
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DO  122  K=1,N2 

IF( I. EQ. K)  THEN 
P2( I,K =  1. 0  - 
ELSE 

P2(  I,K)=0.  0 
END  IF 

122  CONTINUE 
121    CONTINUE 

DO  1121  1=1. N2 

DO  1122  K=1,N2 

IF( I. EQ. K)  THEN 
P3( I,K)=CON2 
ELSE 

P2( I,K)=0. 0 
END  IF 
1122       CONTINUE 
1121    CONTINUE 
C      CALL  FAMILY(  N2 , N2 ,  ' P2 '  ) 
C      CALL  RESULT(N2,N2, I.K.P2) 

CALL  COPYM(N2/N2/ I/L/AB/P2) 
DO  123  1=1 ,N2 
Ul(  I,M)=0.  0 

123  CONTINUE 

CALL  TRANS(N2,M,I,K,F1,FT) 
WRITE( 6.702) 

702  FORMATS  T. 'ENTER  Q(Z)  IN  DESCENDING  ORDER  OF  Z') 
WRITE(6.70$) 

703  FORMATf   T/HIGHEST  COEFF.  MUST  BE  1.  '  ) 
WRITE(  6.704) 

704  FORMATS  '  T ,  f *REMARK:  Q(  Z )  SHOULD  BE  STABLE  POLYNOMIAL') 

DO  157  t=l,N 

WRITE( 6.705 )N-I 

705  FORMAT (  '  T,fCOEFF.  OF  Z**(  ' , 12 ,  '  )  =  '  ) 
READ(  5>)  C3(  I) 

157       CONTINUE 

WRITE( 6,706) 

r,f ENTER  PSTAR(Z)  IN  ASCENDING  ORDER 
OF  Z'  ) 
707) 
T, 'highest  COEFF.  MUST  BE  1.  '  ) 
,706) 
708       FORMAT (  '  T ,  ' * REMARK- 1:  PSTAR(  Z)  IS  THE  DENOMINATOR 

POLYNOMIAL1  ) 
709) 

T. 'OF  MODEL' ) 
7l6) 
710       FORMAT ( '  T, ' *REMARK-2: PSTAR( Z)  SHOULD  BE  STABLE 

POLYNOMIAL1  ) 
DO  158  1=1, N 

WRITE( 6/711)N-I 


706  FORMAT( 

WRITE(6 

707  FORMATf 
WRITE( 6 


WRITE( 6 
709       FORMAT 


WRITE( 


\ 
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711  FORMAT ( '  '  'COEFF.  OF  Z**( ',12, ')=') 

READ(5,*)  C4( I) 
158       CONTINUE 

DO  241  1=1.20 

V(  n=vp(i) 

241       CONTINUE 

DO  114  J=N1,L8 

CALL  CALCUL( N2 , M, N2 . I , L , K, VI , P2 , Fl ) 

CALL  CALCUL(M,M,N2, I ,  L,  K,  V2  ,  FT,  VI ) 

TSP=C0N1+V2(M,M) 

CALL  CALCUL(M,M,N2, I ,L,K, V4,FT,U1) 

TSF=Y( J)-V4  M,M 

TSN=TSF/TSP 

CALL  CALCUL(N2,M,N2. I , L, K, P4, Fl , P3 ) 

CALL  SCALAR  N2,M. I, K,TSN,P4, V5 ) 

CALL  SUM(N2,  M,  I  ,Ky  Ul  ,V5,U2) 

DO  117  1=1, N2 
Ul(  I ,M)=0.  0 
117       CONTINUE 


DO  148  1  =  1, N 
U( I)=V( t) 
148       CONTINUE 

DO  211  1=1,  N 
D( I)=U1(  I,M) 

211  CONTINUE 

DO  212  1  =  1, N 
L4=N+ I 
B(  IJ=U1(L4,M) 

212  CONTINUE 

DO  7011  1  =  1, N 
M9=N-I+1 
Ql( I)=C3(M9) 
7011      CONTINUE 

01(N1)=1.0 

Q2(N1]=0.  0 

DO  4001  1=1, N 

Q3( I)=D( I)-C4( I) 
4001      CONTINUE 

DO  7007  1=1, N 
M8=N-I+1 
Q2( I)=Q3(M8) 
7007      CONTINUE 

DEG1=N+1 

DEG2=2*N 

DEG3=DEG2+2 

DO  6001  I=1,DEG3 
R( I)=0. 0 
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DO  6002  K=1,DEG2 

IF(  (  (I-K).  LT.  0.  0).  OR.  (  (  I-K).  GT.DEG1)  ) 

GO  TO  6002 
IF((  I.  GT.  1).  OR.  (K.  GT.  1})  GO  TO  6503 
6503  R( I-1)=R(  I-1)+Q1(K-1)*Q2(  I-(K-l) ) 

6002         CONTINUE 
6001      CONTINUE 

K9=N6-I+1 

F(IJ=R(K9) 
9003      CONTINUE 
N4=2*N+1 
M1=N+1 

DO  41  J4=1/M1 
DO  51  1=1, N4 

K=I-J4 

A( I, J4i=0. 

IF  (K.  EQ.  0.  )  THEN 

END  IF 

IF  ((K.  GT.  0.  ).  AND.  (  K.  LT.  Ml))  THEN 


A(I, J4)=D(K) 
)  I" 


END  IF 
51        CONTINUE 
41     CONTINUE 

DO  61  J4=1,N 

DO  71  1=1, N4 

A( I,N+1+J4)=0. 

K=I-J4-1 

IF  ((K.  GT.  0.  ).  AND.  (  K.  LT.  Ml)  )  THEN 

A(I/N+i+j4)=b(k) 

END  IF 
71        CONTINUE 
61     CONTINUE 

DO  10  C0LUM=1,N4 

SUM2=0. 

DO  70  J4=COLUM,N4 

SUM2=SUM2+( A( J4,C0LUM) )**2 
70     CONTINUE 

SIGMA=SQRT(SUM2) 

ASIGMA=ABS(  SIGMA) 
IF  (ASIGMA. LT. 0.  00001)  THEN 
WRITE(_6.714) 
714  FORMAT ( i     T, 'MATRIX  A  IS  SINGULAR1) 

END  IF 
DO  80  J4=C0LUM,N4 
G(  J4J=0. 
80     CONTINUE 

G(COLUM)=SIGMA 


113 


SUM2=0. 

DO  90  J4=C0LUM,N4 

SUM2=SUM2+( A( J4, COLUM) -G( J4) ) **2 
90     CONTINUE 

DO  100  J4=C0LUM,N4 
DO  110  K=C0LUM,N4 

IF  ( ABS( SUM2). LT. 0. 00001)  THEN 

TEMP=0. 
ELSE 

TEMP=-2*( A( J4, COLUM) -G( J4) ) *( A( K, COLUM) 
+      ««  TW   "G(k))/SUM2 
END  IF 
IF  ( J4.  EQ. K)  THEN 

TEMP=I. +TEMP 
END  IF 

H( J4,K)=TEMP 
110    CONTINUE 
100    CONTINUE 

DO  131  J4=C0LUMV N4 
DO  130  K=COLUM,N4 
TEMP=0. 
DO  140  L=C0LUM,N4 

TEMP=TEMP+H( J4 , L ) * A(  L , K ) 
140         CONTINUE 

CrJ4,K)=TEMP 

130  CONTINUE 

131  CONTINUE 

DO  150  J4=COLUM7N4 
TEMP=0. 
DO  160  L=COLUM,N4 

TEMP=TEMP+H( J4 , L ) *F( L ) 
160      CONTINUE 

T{  J4 r)  =TEMP 
150    CONTINUE 

DO  170  J4=C0LUM,N4 
DO  180  K=C0LUM,N4 
A(  J4,K)=C( J4/K) 
180      CONTINUE 

F( J4)=T( J4) 
170    CONTINUE 
10     CONTINUE 

X(  N4 )  =F(  N4 )  /A(  N4 ,  N4 ) 
I=N4-1 
812      SUM2=0. 
M2=I+1 
DO  30  J4=M2.N4 

SUM2=SUM2  +  A(  I, J4)*X( J4) 
30       CONTINUE 

XLI)=LF( I)-SUM2)/A( 1,1) 

IF  (ABS(X(I)).LT.  0.00001)  THEN 
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X(  I  =0. 
END  IF 
1  =  1-1 

IF( I. GE. 1)  GO  TO  812 
DO  141  1=1. Nl 

S2( I)=X(I) 

141  CONTINUE 

DO  142  I=N3,N4 
N5=I-N1 
S1(N5)=X(  I) 

142  CONTINUE 
UP(  J)=0.  0 

DO  5001  1  =  1, N 

UP("ji=UP^JJHS2(^L7j-C3(I))*U(J-I)  +  Sl(I)*Y(J-I) 

5001     CONTINUE 

U(J)=(UP(J)+V(J))/(1-S2(1)) 
J5=J+1 
Y(J5)=0.  0 
DO  3003  1  =  1, N 
M6=N1-I 
Y( J5)=Y( J5)-C6( I)*Y( J5-M6)+C5( I ) *U( J5-M6) 

3003  CONTINUE 

DO  189  1=1, Nl 

YM(  1 1=0.  0 
189     CONTINUE 

YM(  J5]=0.  0 
DO  2508  I  =  1,N 

M8=N-I+1 

Bl(  I)=C5(M8) 
2508    CONTINUE 

DO  3004  1=1, N 

YM( J5)=YM( J5)-C4( I )*YM( J5-I )+Bl(  I)*V(  J5-I) 

3004  CONTINUE 

MK( J5]=-Y( J5) 
DO  3005  1=1. N 

MK(  J5)=MK( J5)+B(  I ) *V( J5-I ) -C4(  I)*Y( J5-I) 

3005  CONTINUE 

MU(  J5)=ABS(MKCJ5) ) 
IF(MU(J5J.  LT.  20.  0   THEN 

V(J5)=fc(J5) 
ELSE 

V( J5)=VC( J5)+VP( J5) 
END  IF 

DO  120  I=1,N2 
Fl(  I,M)=0.  0 
120       CONTINUE 

DO  118  I=1,N 
L1=N+I 
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118 
119 

3001 


3002 

114 

2255 


8883 
8882 


1111 

C 

8002 
976 


1025 


J2=J+1-I 

Fl( I.M)=-Y(  J2) 

Fl(Ll,M)=U(  J2) 
CONTINUE 
DO  119  1=1, N2 

FT(M,  I)=0.  0 
CONTINUE 


DO  3001  1=1. N2 

U4( I,M)=6. 0 
CONTINUE 

CALL  SUB(N2/N2/I/L/U3/U1/U4) 
L6=J-N 
C7(L6]=0.  0 
DO  3002  I=1,N2 

C7(L6)=C7(L6)+U4(  I,M)*U4(  I,M) 
CONTINUE 

C8(  L6)=SQRT( C7( L6) ) 
CONTINUE 

WR I TE (8 , 2  2  5  5 ) ' PARAMETER  ERROR  * 
FORMATS 12X.A15) 
DO  8882  1=1,100 

WRITE(8,8883)  C8( I) 

FORMAT (  15X,F12. 6) 
CONTINUE 


WRITE( 8.1111) 'SYS.  OUTPUT' 

F0RMAT( 9X, Al 1 , 6X , Al 1 ) 

DO  976  1=1,100 

WRITE  (6,8002)  Y(I),YM(I 

WRITE  (8.8002)  Y(I),YM(I 

FORMAT? 9X,F 10.  6,6X,F10.  6 

CONTINUE 

DO  1025  I=l/L8 

SE( i;2]=YM( I) 
KD(  I >=I 

CONTINUE 

65827 

65827 

65827 

65827 

65827 
TL0T(KD,SE,L8.2, 'TIME 
V1PL0T(  KD,C8,L8,  1,  TTIME 
V2PL0T  KD,MU, L8. 1, 'TIME 
V3PL0T  KD, 7,18, 1, 'TIME 
V4PL0T  KD,U,L8, 1,  'TIME 


' MOD.  OUTPUT 
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STOP 
END 
C   *THIS  SUBROUTINE  PRINTS  THE  ARRAYS  AS  AN  MATRIX  FORM. * 
SUBROUTINE  RESULT( H, 0, P, S, T) 
INTEGER  H  ,  0  ,  P  ,  S 
REAL  T  (  H  ,  0  ) 
DO  80   P  =  1  ,  H 

WRITE( 8.70)(T(P,S) ,S=1,0) 
C         PRINT  70.  (T(_P,SJ,S=1.0) 
70        FORMAT(  T0\T2,9X,10(3X,F12.  6)  ) 
80      CONTINUE 
RETURN 
END 
C   ***THIS  SUBROUTINE  CALCULATES  THE  MULTIPLICATION  OF 
TWO  MATRICES. *** 
SUBROUTINE  CALCUL( U.  V, Y, Z, X, ZX,  ZT,  W,  Q) 
INTEGER  U,V,Y, Z.X, ZX 
REAL  ZT(y,v},W(U,Y),Q(Y,V) 
DO  93  Z=1,U 
DO  92  X=1,V 

ZT(Z.X)  =  0.  0 
DO  91  ZX=1,Y 

ZT( Z,X)=ZT(  Z,X)+W( Z,ZX)*Q( ZX,X) 

91  CONTINUE 

92  CONTINUE 

93  CONTINUE 
RETURN 
END 

c  ***THis  SUBROUTINE  WRITES  THE  MATRIX  NAME  AND  ROW, 
COLUMN  NUMBERS. ***/ 

SUBROUTINE  FAM I LY( ROW, COLUMN, NAME) 

INTEGER  ROW, COLUMN 

CHARACTER* 2  NAME 

WRITE( 8,94) 
C       PRINT  94,  rMi 
94'     FORMATfTl ,  *  0 '  T2 , 10X, T12 . A6 , T19 , IX, T20 , A2 ) 

WRITE(8,95)  ' ROW  NUMBER.  f , ROW 
C       PRINT  95.TROW  NUMBER.  T, ROW 

95  FORMATfTl/O'  ,T2,  10X,  T12  ,  All,  T23 ,  8X,  T31,  12) 
WRITE(8,96)  'COLUMN  NUMBER.  .COLUMN 

C       PRINT  96,rCOLUMN  NUMBER. T. COLUMN 

96  FORMAT( Tl , r0r , T2 , 10X, T12 , Al4, T26 , 5X, T3 1 , 12 ) 
RETURN 
END 

c   ***THIS  SUBROUTINE  CALCULATES  THE  SUMMATION  OF  THE 
C      absolute  values  OF  THE  SAME  COLUMN  ELEMENTS. *** 

SUBROUTINE  ALI( I 1 , 12 , HI , H2 , Gl , W2 ) 

INTEGER  II, 12, HI. H2 

REAL  W2( 1, I2),Gl( II. 12) 
DO  41  H2=l, 12 


,94)  ^MATRIX*  NAME 

4. rMATRIX^ .NAME 

Tl,  *  0T,T2 , 10X, T12. A6, 

nt\       i  r>/""\r„7    MTTiurotrr)  nr\ 
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DO  42  Hl=l. II 
W2(  1,H2)=W2(  1/H2)+ABS(G1(H1,H2) ) 

42  CONTINUE 
41      CONTINUE 

RETURN 
END 
C   ***THIS  SUBROUTINE  MULTIPLIES  THE  MATRIX  BY  SCALAR 
NUMBER. *** 
SUBROUT INE  SCALAR(  Kl , K2 , Gil , G2 , PI , G3 , G4) 
INTEGER  K1,K2,G11,G2 
REAL  G3(K1,K2KG4(K1/K2)/P1 
DO  43  G2=1,K2 

DO  44  G11=1,K1 
G4(G11,G2)=G3(G11,G2)*P1 

44  CONTINUE 

43  CONTINUE 
RETURN 
END 

C   ***THIS  SUBROUTINE  CALCULATES  THE  SUM  OF  THE  TWO 
MATRICES. *** 
SUBROUTINE  SUM(K3.K4, 15, K5,X1,X2 ,X3 ) 
INTEGER  K3,K4,I5.K5 

REAL  XI ( K3 , K4J , X2 ( K3 , K4 ) , X3  (  K3 , K4 ) 
DO  45  ^5=1^4 

DO  46  15=1. K3 

X3(  I5/K5)=X1( I57K5)+X2( I5,K5) 
46      CONTINUE 

45  CONTINUE 
RETURN 
END 

C   ***THIS  SUBROUTINE  COPPIES  THE  MATRICES.*** 
SUBROUTINE  COPYM(K3/K4/ I5,K5/X1/X2) 
INTEGER  K3,K4, I5.K5 
REAL  X1(K3/K4),X2(K3/K4) 
DO  45  K5=1#K4  . 

DO  46  15=1. K3 

Xl( I5#K5)=X2( I5,K5) 

46  CONTINUE 
45      CONTINUE 

RETURN 
END 
C   *THIS  SUBROUTINE  FINDS   THE  TRANSPOSE   OF  THE  MATRIX.* 
SUBROUTINE  TRANS(K3,K4/ I5/K5,X1,X2) 
INTEGER  K3,K4, I5.K5 
REAL  X1£K3 , K4J , X2 ( K4 , K3 ) 


46      CONTINUE 
45      CONTINUE 


DO  45  K5=1,K4 

DO  46  15=1, K3 
X2(K5, I5)=X1( I5,K5) 
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RETURN 
END 
C   ***THIS  SUBROUTINE  CALCULATES  THE  SUB.  OF  THE 
TWO  MATRICES.  *** 
SUBROUTINE  SUB(K3.K4, 15 , K5 , XI,  X2 , X3 ) 
INTEGER  K3,K4,I5,K5 

REAL  Xl(  K3 , K4J ,X2( K3 , K4) , X3(  K3 , K4) 
DO  45  K5=1,K4 

DO  46  15=1. K3 

X3( I5,K5)=X1( I5/K5)-X2(  15, K5) 
46      CONTINUE 
45      CONTINUE 
RETURN 
END 

C\         •k-k-k-kk-k-kkkkk'kk-kk'kkk-kkk'kk'k'kkickkjck'k'kit'kick'kk-k'kk 

C   *  THIS   SUB.  PLOTS   THE   ARRAYS   USING    * 
C   *  THE  DISSPLA  PROGRAM.  * 

rj      ******************************************** 

C 

SUBROUTINE  VPLOT(  T,  U,M,N,  LABELX,  LABELY) 

H-kk-kkkkkkkkkk-k 

REAL  U( 800.4) ,T( 800) ,Y( 800) .YX(4) ,YN( 4) 
INTEGER   PMAX, PMIN, M, N, LABELX, LABELY 
C 

CALL  AMAX(T,M,TMAX,PMAX) 
CALL  AMINfT^TMIN/PMIN) 
DO   20   J   =  1,N 

DO  10  I  =  1,M 

Y(I)  =U(I,J) 
10       CONTINUE 

CALL  AMAX( Y . M , YMAX , PMAX ) 

YX(  J  )    =  YmAX 

CALL  AMIN(Y,M,YMIN,PMIN) 

YN(J)  =YMIN 

WRITE(6,11)  YMAX,YMIN 
20        CONTINUE 

WRITE(6,13)(YX(  I), 1  =  1 ,4) 
WRITE(  6  .  13 j  (  YN(  ij  .  1  =  1.  4) 
CALL  AMAX(YX,N, YMAX, PMAX) 
CALL  AMIN(YN,N,YMIN,PMIN) 
WRITE(6. 11)  YMAX,YMIN 
CALL  TEK  618 
CALL  BLOWUPfl. 2) 
CALL  PAGE( 11.  ,8.  5) 
C      CALL  NOBRDR 

CALL  AREA2D(9.  ,6.  ) 

CALL  XNAME( LABELX, 6) 

CALL  YNAME(  LABELY, 6) 

CALL  HEADIN  ( TPLANT  AND  MODEL  OUTPUTS  $' ,100,1. ,1) 

CALL  CROSS 


119 


CALL  GRAF(TMIN, 'SCALE' ,TMAX,YMIN, 'SCALE' ,YMAX) 
CALL  SPLINE 
DO  40  J  =  1,N 

DO  30  I  =1,M 

Y(I)  =u(l,J) 
30       CONTINUE 

IF(J.  EQ.  2)  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 
40    CONTINUE 

CALL  ENDPL(O) 
C      CALL  DONEPL 

FORMAT^  \10X,2G12.  5/) 
FORMAT)'  ' ,10X,4G12. 5/J 
RETURN 
END 
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SUBROUTINE  V1PL0T(  T,U,M,  N,  LABELX,  LABELY) 

C************* 

REAL  U( 800.4) ,T( 800) ,Y( 800) , YX( 4 ) ,YN( 4) 
INTEGER  PMAX, PMIN, M, N, LABELX, LABELY 
C 

CALL  AMAX(T,M,TMAX,PMAX) 
CALL  AMIN(T,M,TMIN,PMIN) 
DO  20  J  =  1,N 

DO  10  I  =1.M 

Y(I)  =u(l,J) 
10       CONTINUE 

CALL  AMAX(  Y , M , YMAX , PMAX ) 

YX(J)  =  YMAX 

CALL  AMIN(Y,M,YMIN,PMIN) 

YN( J)  =YMIN 

WRITE(6,11)  YMAX,YMIN 
20    CONTINUE 

WRITE( 6,13)(YX( I), 1=1,4) 
WRITECe.lSJlYN  ij  . 1=1 ,  4J 
CALL  AMAX(YX,N, YMAX, PMAX) 
CALL  AMIN(YN,N,YMIN,PMIN) 
WRITE(6.  11)  YMAX,YMIN 
CALL  TEK  618 
CALL  BL0WUP(1. 2) 
CALL  PAGE( 11. ,8. 5) 
C      CALL  NOBRDR 

CALL  AREA2D(9.  ,  6.  ) 

CALL  XNAME( LABELX, 6 ) 

CALL  YNAME( LABELY. 6J 

CALL  HEADIN  ( TPARAMETER  ERROR  $' ,100,1. ,1) 

CALL  CROSS 

CALL  GRAF(TMIN, "SCALE' ,TMAX,YMIN, 'SCALE' ,YMAX) 

CALL  SPLINE 

DO  40  J  =1,N 


120 


DO  30  I  =1.M 

Y(I)  =U(I,J) 
30       CONTINUE 

IF(J.  EQ.  2)  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 
40    CONTINUE 

CALL  ENDPL(O) 
C      CALL  DONEPL 

FORMAT(  '  7 , 10X,2G12.  5/) 
FORMAT  '  ' , 10X, 4G12.  5/j 
RETURN 
END 


11 
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SUBROUTINE  V2PL0T(  T,U,M,N, LABELX, LABELY) 

REAL  U( 800.4) ,T( 800) , Y( 800 ) . YX( 4) . YN(  4) 
INTEGER  PMAX, PMIN, M, N, LABELX, LABELY 
C 

CALL  AMAX(T,M,TMAX,PMAX) 
CALL  AMIN[T,M/TMIN/PMIN) 
DO  20  J  =1,N 

DO  10  I  =1.M 

Y(I)  =U(I,J) 
10       CONTINUE 

CALL  AMAX( Y . M , YMAX , PMAX ) 

YX(  J )    =  YMAX 

CALL  AMIN( Y,M,YMIN,PMIN) 

YN(  J)  =YMIN 

WRITE(6,11)  YMAX,YMIN 
20    CONTINUE 

WRITE( 6,13)(  YX(  I), 1  =  1,4) 
WRITE?  6, 13j( YN(  ij , 1  =  1, 4J 
CALL  AMAX(  YX, N, YMAX, PMAX) 
CALL  AMIN(YN,N.YMIN,PMIN) 
WRITE(6. 11)  YMAX,YMIN 
CALL  TEK  618 
CALL  BL0WUP(1.  2) 
CALL  PAGE( 11.  ,8.  5) 
C      CALL  NOBRDR 

CALL  AREA2D(9.  ,6.  ) 

CALL  XNAME( LABELX, 6) 

CALL  YNAME  LABELY, 6j 

CALL  HEADIN  ( TOUTPUT  PREDICTION  ERROR  $',100,1. ,1) 

CALL  CROSS 

CALL  GRAF(TMIN, 'SCALE' ,TMAX,YMIN, "SCALE' ,YMAX) 

CALL  SPLINE 

DO  40  J  =1,N 

DO  30  I  =1.M 

Y(IJ  =U(I,J) 
30       CONTINUE 
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IF(J.  EQ.  2)  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 
40    CONTINUE 

CALL  ENDPL(O) 
C      CALL  DONEPL 

11   FORMAT^   ,10X,2G12.5/ 
13   FORMAT (  *  ' ,10X,4G12.5/ 
RETURN 
END 
C 

SUBROUTINE  V3PLOT( T,U, M,N, LABELX, LABELY) 

REAL  U(  800.4) ,T( 800) , Yf 800) . YX(4)  ,YN( 4) 
INTEGER  PMAX, PMIN, M, N, LABELX, LABELY 
C 

CALL  AMAX(T,M,TMAX,PMAX) 
CALL  AMIN(T,M,TMIN,PMIN) 
DO  20  J  =  1,N- 
DO  10  I  =1.M 

Y(I)  =u(l,J) 
10       CONTINUE 

CALL  AMAX(Y.M,YMAX,PMAX) 
YX( J )  =  YMAX 
CALL  AMIN(  Y/M/YMIN/PMIN) 
YN(  J)  =YMIN 
WRITE(6,11)  YMAX,YMIN 
20    CONTINUE 

WRITE( 6/13)(YX( I), 1=1,4) 
WRITE(  6.  13  j  (  YN(  IJ  ,  1  =  1,  4J 
CALL  AMAX(YX,N, YMAX, PMAX) 
CALL  AMIN(YN.N,YMIN,PMIN 
WRITE(_6.  11]  YMAX,YMIN 
CALL  TEK  618 
CALL  BLOWUP(l.  2) 
CALL  PAGE( 11.  ,8.  5) 
C      CALL  NOBRDR 

CALL  AREA2D(9.  ,6.  ) 

CALL  XNAME( LABELX, 6 ) 

CALL  YNAME  LABELY, 6 j 

CALL  HEADIN  ( TEXTERNAL  INPUT  $',100,1.  ,1) 

CALL  CROSS 

CALL  GRAF(TMIN, 'SCALE' ,TMAX,YMIN, 'SCALE' ,YMAX) 

CALL  SPLINE 

DO  40  J  =1,N 

DO  30  I  =1,M 

Y(I)  =U(I,J) 
30       CONTINUE 

IF(J.  EQ.2)  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 
40    CONTINUE 
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CALL  ENDPL(O) 
C      CALL  DONEPL 

11   FORMAT^  7  ,10X,2G12.  5/) 
13   FORMAT (  '  '  /10X/4G12.  5/j 

RETURN 

END 

SUBROUTINE  V4PLOT( T,U,M, N, LABELX, LABELY) 

REAL  U(800.4),T(800).Y(800),YX(4),YN(4) 

INTEGER  PMAX, PMIN, M, N,  LABELX, LABELY 
C 

CALL  AMAX(T/M/TMAX/PMAX) 
CALL  AMIN(T,M/TMIN/PMIN) 

DO  20  J  =  1,N 

DO  10  I  =1,M 

Y(I)  =U(I,J) 
10       CONTINUE 

CALL  AMAX(Y.M,YMAX,PMAX) 
YX  (  J )    =  YM  AX 
CALL  AMINf  Y,M,YMIN,PMIN) 
YN( J)  =YMIN 
WRITE(6,11)  YMAX,YMIN 
20    CONTINUE 

WRITE( 6,13)(YX( I), 1=1,4) 
WRITE(  6,  13j(  YN(  ij  ,  1  =  1,  4) 
CALL  AMAX(YX,N,YMAX,PMAX) 
CALL  AMIN(YN,N.YMIN,PMIN) 
WRITE( 6^111  YMAX,YMIN 
CALL  TEfc  618 
CALL  BLOWUP(l.  2) 
CALL  PAGE( 11.  ,8.  5) 
C      CALL  NOBRDR 

CALL  AREA2D(9.  ,6.  ) 

CALL  XNAME( LABELX, 6 ) 

CALL  YNAME( LABELY, 6 ) 

CALL  HEADIN  (T INPUT  TO  THE  PLANT  $',100,1. ,1) 

CALL  CROSS 

CALL  GRAF(TMIN, 'SCALE' ,TMAX,YMIN, 'SCALE' ,YMAX) 

CALL  SPLINE 

DO  40  J  =1,N 

DO  30  I  =1,M 

Y(IJ  =U(I,J) 
30       CONTINUE 

IF(J.  EQ.2)  CALL  CHNDOT 
CALL  CURVE  (T,Y,M,0) 
40    CONTINUE 

CALL  ENDPL(O) 
C      CALL  DONEPL 

11   FORMAT^  7 ,10X,2G12.  5/) 
13   FORMAT}'  ' ,10X,4G12. 5/j 
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RETURN 

END 
C 

SUBROUTINE  AMAX(  Y, N, YMAX, PMAX) 
C 

Q  ******************************************** 

C         *   SUBROUTINE  TO  COMPUTE  THE  LARGE  ELELEMENT 
C         *   IN  A  GIVEN  ROW  OR  COLUMN  IN  ARRAY  TAT. 

C\  ******************************************** 

c 

C         ***  VARIABLE  DECLARATION  *** 
C         REAL  AMAX 

INTEGER  PMAX 
DIMENSION  Y(N) 
C 

YMAX=Y(1) 

DO  5100  1=2, N 

IF(Y(  I).  LE.  YMAX)  GO  TO  5000 
YMAX=Y(  I ) 
PMAX=I 
5000         CONTINUE 
5100      CONTINUE 
C         MAXROW=POS 
RETURN 
END 
C 
C 

SUBROUTINE  AMIN(  Y,N,  YMIN,  PMIN) 
C 

Q  ******************************************** 

C         *   SUBROUTINE  TO  COMPUTE  THE  SMALL  ELEMENT 

C         *   IN  A  GIVEN  ROW  OR  COLUMN  IN  ARRAY  TAT. 

rj         ******************************************** 

C 

c  -       ***  VARIABLE  DECLARATION  *** 

C         REAL  AM IN 

INTEGER  PMIN 
DIMENSION  Y(N) 
C 

YMIN=YC1) 

DO  6100  1=2 ,N 

IF( Yfl). GE.  YMIN)  GO  TO  6000 
YMIN=Y(  I ) 
PMIN=I 
6000  CONTINUE 

6100      CONTINUE 
RETURN 
END 
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